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RECOVERY  OF  FROZEN  FLOW  LOSSES  IN  ARCJETS 


Grant  No.  AFOSR-91-0318 

Final  Technical  Report 
(1991-1993) 

V.  V.  Subramaniam,  V.  Babu,  and  S.  Aithal 
Department  of  Mechanical  Engineering 
The  Ohio  State  University 
Columbus,  Ohio  43210 


1.0  INTRODUCTION 

The  recent  successful  commercial  launch  of  two  1.8  KW  hydrazine  arcjet  thrusters  on 
AT&T’s  TELSTAR  IV  satellite  marks  a  new  era  in  electric  propulsion.  The  arcjet’s  key 
contribution  in  the  commercial  sector  is  in  reduced  launch  costs  and  extended  satellite  operational 
lifetime.  This  is  achieved  primarily  by  increasing  the  specific  impulse  above  that  of  conventional 
chemical  thrusters.  For  military  applications  however,  launch  costs  are  a  secondary  concern 
compared  to  extended  satellite  lifetime  for  multiple  missions  and  quick  response  to  mission  goal 
changes.  These  military  applications  will  require  arcjets  capable  of  operation  over  a  wide  range  of 
powers  from  ~1  KW  up  to  ~50  KW.  Present  understanding  of  arcjet  physics,  arcjet  design,  and 
arcjet  performance  especially  with  regard  to  scale-up  with  input  power  level.  This  report  describes 
progress  made  during  the  first  two  years  of  a  four-year  research  program  designed  to  address  these 
issues. 

The  arcjet  thruster  is  one  in  a  class  of  electrothermal  devices  that  generates  thrust  by  ohmic 
heating  of  the  subsonic  propellant  stream,  followed  by  rapid  expansion  to  supersonic  speeds.  An 
important  distinction  between  these  thrusters  and  conventional  chemical  thrusters  is  the  fact  that  the 
input  electrical  energy  is  transferred  to  the  propellant  via  Ohmic  or  Joule  heating,  as  opposed  to 
combustion.  However,  the  arcjet  shares  a  key  feature  in  common  with  its  chemical  thruster 
counterparts,  namely,  frozen  flow  losses. 

Frozen  flow  losses  refers  to  the  fraction  of  input  energy  that  is  added  to  the  propellant  gas 
flow  which  does  not  appear  in  the  form  of  exit  kinetic  energy  (i.e.  translational  energy).  The 
thrust  generated  derives  solely  from  this  exit  kinetic  energy.  Fig.  1  shows  how  the  input  electrical 
energy  is  partitioned  into  several  channels.  It  is  clear  that  only  a  fraction  of  the  input  energy 
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surfaces  as  useful  work  done  by  the  thrust  force.  Frozen  flow  losses  include  penalties  associated 
with  dissociation  of  propellant  molecules,  ionization  processes,  electronic  excitation  of 
molecules/atoms,  and  vibrational/rotational  excitation  of  molecules.  Each  of  these  processes  has  a 
characteristic  relaxation  time  associated  with  it,  and  the  resulting  non-equilibrium  effects  are 
drastically  coupled  to  the  flow  field.  The  term  “non-equilibrium”  has  been  freely  given  different 
meanings  in  the  electric  propulsion  community.  For  this  reason,  the  type  of  non-equilibrium 
phenomena  of  relevance  to  arcjet  flows  is  clarified  in  the  following  section,  followed  next  by  a 
description  of  completed  and  on-going  research. 

This  final  report  summarizes  the  research  performed  under  grant  no.  91-0318,  during  the 
period  September  1991  -  September  1993. 

2.0  BACKGROUND 

Although  early  tests  of  arcjet  thrusters  in  the  1960s  contributed  to  the  development  of  a 
prototype  device[l-4],  dangerous  halting  of  further  development  in  the  1970s  with  subsequent 
revival  of  space-based  defense  concepts  in  the  1980s  have  led  to  stagnation  in  fundamental 
understanding  of  the  operating  characteristics  and  performance  of  these  devices.  In  the  mid-1980s 
therefore,  the  1  KW  and  30  KW  designs  from  the  1960s  were  revived,  tested  for  lifetime  and 
endurance  using  ammonia  as  the  propellant,  and  appeared  to  meet  the  mission  requirement 
(ISp~800s)  set  by  the  SDIO[5].  However,  in  1990,  systematic  errors  were  discovered  in  the  thrust 
measurements  made  at  JPL  which  led  to  performance  evaluations  that  failed  to  meet  the  mission 
criteria  for  the  ammonia  arcjet.  Consequently,  all  research  efforts  in  the  U.S.  turned  toward  the 
hydrogen  arcjet  which  could  easily  generate  specific  impulses  on  the  order  of  1500s.  This  work 
was  proposed  at  that  time  with  the  following  simple  argument  in  mind.  Although 
the  ammonia  arcjet  had  an  Isp  far  lower  than  the  hydrogen  arcjet,  propellant  storage  and  storage 
cost  considerations  favor  the  ammonia  arcjet  over  the  hydrogen  arcjet.  Any  gain  in  Isp  in  the 
ammonia  arcjet  would  therefore  be  a  benefit,  in  addition  to  the  ease  of  storage  and  lowering 
propellant  storage  costs.  Unfortunately,  polatomic  propellants  such  as  ammonia  usually  react  to 
form  significant  amounts  of  highly  vibrationally  excited  diatomic  species.  Thus,  a  substantial 
amount  of  the  input  electrical  energy  is  lost  in  frozen  flow  (i.e.  rotational,  vibrational,  and 
electronic  mode  non-equilibrium).  Clearly,  if  these  frozen  flow  losses  are  to  be  reduced  and  if 
sustained  operation  is  to  be  achieved  at  high  power  levels  and  high  specific  impulses,  a 
fundamental  understanding  of  the  transport  and  chemical  processes  is  necessary.  Scale-up  in 
power  is  not  possible  without  consideration  of  the  disequilibria  between  internal  and  translational 
modes  of  molecular  motion.  It  is  therefore  important  to  define  exactly  what  is  meant  by  a  “non¬ 
equilibrium”  process. 

The  words  “non-equilibrium”  are  used  in  the  scientific  community  to  denote  a  process 
whose  characteristic  time  scale  is  comparable  to  (i.e.  neither  orders  of  magnitude  larger  nor  smaller 
than)  the  time  scale  set  by  a  flow.  There  are  in  general  two  time  scales  in  a  real  flowing  fluid.  The 
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first  is  the  convective  time  scale  defined  as  a  macroscopic  characteristic  length  (such  as  constrictor 
diameter,  constrictor  length,  channel  diameter,  etc.)  divided  by  a  characteristic  velocity,  and  the 
second  is  a  diffusive  time  scale  defined  as  the  square  of  a  characteristic  length  divided  by  the 
kinematic  viscosity.  According  to  classical  thermodynamics,  a  system  is  in  equilibrium  if  there 
are  no  gradients  in  any  properties  of  the  fluid.  However,  when  there  is  a  flowing  compressible 
fluid,  there  are  gradients  in  velocity,  density,  etc.  Consequently,  one  of  the  “ideal”  flows  that  is 
commonly  used  for  model  comparisons  is  one  that  is  in  local  thermodynamic  equilibrium.  Such 
an  equilibrium  flow  is  one  in  which  the  time  scales  fen*  all  processes  are  infinitely  smaller  than 
either  the  convective  or  diffusive  time  scale.  A  second  “ideal”  flow  is  “ frozen ”  flow,  which 
refers  to  the  case  where  the  time  scale  for  a  particular  process  is  much  longer  than  the  convective  or 
diffusive  time  scale.  Hence,  the  chemical  composition  is  “frozen”  or  unchanging. 

A  given  process  (i.e.  dissociation,  ionization,  recombination,  reaction,  etc.)  is  out  of 
equilibrium  in  a  flow  if  the  time  scale  for  this  process  is  comparable  to  the  convective  time  scale. 
Hence  it  is  possible  for  flows  to  consist  of  many  processes,  some  of  which  are  in  local 
thermodynamic  equilibrium,  some  that  are  frozen,  and  others  that  are  out  of  thermodynamic 
equilibrium.  The  processes  that  are  out  of  equilibrium  in  this  sense  are  called  “non-equilibrium 
processes”.  Since  classical  thermodynamics  does  not  account  for  the  rate  at  which  a  process  may 
take  place,  there  is  no  recourse  but  to  resort  to  the  area  of  chemical  kinetics.  Chemical  kinetics  is  a 
term  broadly  used  by  many  scientists  to  denote  finite-rate  reactions  (as  opposed  to  infinite  rates 
which  correspond  to  equilibrium  conditions,  and  zero  rates  which  correspond  to  frozen  states). 
Thus  any  reaction  (to  be  distinguished  from  an  elementary  process  or  elementary  reaction  )  is 
described  by  some  rate  at  which  it  proceeds.  In  many  applications,  reactions  are  described  by 
overall  rates  and  an  overall  chemical  equation  despite  the  fact  that  the  actual  process  takes  place  via 
several  elementary  steps. 

In  reality,  chemical  processes  (dissociation,  ionization,  recombination,  reaction,  etc.) 
proceed  at  different  rates  depending  on  the  detailed  initial  state  of  the  reactants.  For  atomic 
reactants  this  means  that  the  reaction  may  proceed  from  a  specific  electronic  configuration  or  state 
of  the  atom.  For  molecular  reactants,  the  reaction  may  proceed  from  specific  vibrational  and 
rotational  levels  in  a  specific  electronic  state.  The  rate  of  such  a  process  or  reaction 
therefore  depends  on  vibrational  quantum  number  and  rotational  quantum  number  for  a  given 
electronic  state.  The  kinetic  rates  for  such  elementary  processes  are  what  are  known  as  “ state- 
specific  kinetic  rates”  .  There  has  been  a  rapid  growth  in  the  generation  of  state-specific  rate 
data  for  many  molecules  over  the  past  decade  or  so,  due  in  part  to  the  sophisticated  diagnostics 
methods  that  have  been  devised.  The  time  is  therefore  ripe  to  make  use  of  such  detailed  and 
accurate  information  where  available,  and  to  combine  them  with  recent  rapid  advances  in 
supercomputing. 

For  high  speed  flows  of  interest  to  rocket  engines  (chemical  or  electrical),  such  state- 
specific  information  is  vital  because  it  determines  transport  properties  (since  they  are  dependent  on 
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the  chemical  composition  of  the  gas  mixture),  which  are  continuously  changing  due  to  the  varying 
flow  field.  However,  also  of  vital  importance  is  the  fact  that  when  molecular  species  are  present  in 
high  speed  flows  (and  this  is  true  of  all  molecular  propellants  in  supersonic  nozzle  flows  whether  it 
is  a  chemical  or  an  electrical  thruster),  a  substantial  amount  of  energy  is  tied  up  in  the  internal 
modes  (vibration,  rotation,  electronic).  This  energy  is  of  no  use  in  generating  thrust 
because  it  is  not  in  translation  (i.e.  it  is  not  directed  kinetic  energy).  Consequently,  energy  tied  up 
in  these  internal  modes  is  “lost”,  and  termed  “frozen  flow  losses”.  The  term  “frozen  flow 
losses”  is  apt  because  the  characteristic  relaxation  time  for  the  vibrational  modes,  or  rotational 
modes,  or  electronic  modes  are  longer  than  the  local  convective  time  scale.  Hence,  the  energy  is 
“frozen”  in  these  modes.  This  is  not  to  be  confused  with  the  “frozen  flow”  concept  addressed 
earlier  which  usually  refers  to  non-varying  chemical  composition. 

With  the  aforementioned  definitions  and  distinctions  in  mind,  the  existing  literature  on 
arejet  flows  can  now  be  reviewed.  Several  research  groups  have  undertaken  study  of  arejet  flow 
dynamics  both  computationally  and  experimentally.  Their  work  will  be  reviewed  here  briefly. 
Schrade  et.  al.  [6]  have  modelled  two-dimensional,  axi-symmetric,  fully  ionized  flow  in  an  MPD 
arejet.  While  this  has  relevance  to  the  hybrid  (electromagnetic/electrothermal)  device,  it  is  not 
directly  applicable  to  the  arejet  which  is  a  purely  electrothermal  device.  The  main  drawback  is  that 
the  detailed  finite  rate  chemistry  which  is  vital  to  arejet  performance,  is  not  modelled.  More 
recently,  Schrade  et.  al.[7,8]  have  extended  their  two-fluid,  quasi  one-dimensional  model  in  order 
to  study  the  interaction  between  the  flow  in  the  hot  arc  and  the  cooler  outer  flow  in  the  constrictor 
region  of  the  arejet.  Their  focus  has  been  mainly  to  understand  the  arc  attachments  and  arc 
dynamics  in  the  constrictor  region[9].  Butler  et.  al.[10].  King  and  Butler[ll],  and  Rhodes  and 
Keefer[12]  have  also  begun  modelling  of  two-temperature  axi-symmetric  flows  in  arejets.  A 
common  denominator  however  in  all  the  aforementioned  works  is  that  the  state-resolved  chemistry 
(especially  including  vibrational  and  electronic  non-equilibrium)  is  neglected.  Additionally,  many 
of  these  investigators  are  utilizing  methods  (such  as  SIMPLE[13],  for  instance)  commonly  avoided 
for  viscous  supersonic  internal  flows,  on  grids  that  are  coarse  (less  than  50  x  50  for  instance). 
Rhodes  and  Keefer  even  mention  that  when  swirl  is  added  to  their  model,  no  converged  solution 
could  be  obtained[12].  Such  results,  although  first  in  modelling  of  arejet  flows,  should  therefore 
be  treated  with  cautious  optimism. 

More  recently,  several  works[14-16]  have  attempted  to  improve  the  state-of-the-art  in 
modelling  of  arejet  flows.  However,  these  efforts  have  been  plagued  by  problems  ranging  from 
long  run  times  (on  the  order  of  hundreds  of  CPU  hours  on  CRAY  supercomputers)  to  ad-hoc 
assumptions  such  as  the  use  of  conductivity  floors.  Run  times  on  the  order  of  hundreds  of  CRAY 
hours  hardly  constitutes  a  useful  design  tool,  and  ad-hoc  fixes  such  as  the  use  of  a  “conductivity 
floor”  will  not  be  very  useful  in  yielding  predictive  design  tools  (since  adjustable  constants  can  be 
made  to  vary  to  fit  almost  any  experimental  result). 

In  contrast,  the  method  developed  in  the  work  reported  here  is  found  to  be  highly  accurate 
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and  reliable.  Furthermore,  our  scheme  is  equally  applicable  to  quasi  1-D,  2-D  axi- symmetric,  as 
well  as  3-D  problems  including  large  sets  of  master  rate  equations  describing  rate  processes.  The 
ultimate  power  of  this  technique  is  evident  especially  when  applied  on  modem  supercomputers. 
The  present  research  therefore  represents  a  significant  advance  in  the  state-of-the- 
art. 

3.0  RESEARCH  OBJECTIVES 

The  objectives  of  this  four-year  research  program  are  (1)  to  study  and  quantify  detailed 
vibrational  and  electronic  non-equilibrium  energy  transfer  processes  in  the  arcjet  nozzle,  (2)  to 
quantify  the  ro-vibrational  populations  using  modem  information  on  state-resolved  rate  constants 
to  enable  accurate  determination  of  frozen  flow  losses,  and  (3)  numericallly  study  the  effects  of 
dilution  of  the  propellant  feed  with  a  fast  Vibration-Translation  (VT)  relaxer,  on  recovery  of  frozen 
flow  losses. 

4.0  STATUS  OF  RESEARCH  EFFORT 

The  discussion  below  summarizes  the  highlights  of  the  research  efforts  under  the  first  year 
of  this  grant.  A  more  detailed  discussion  can  be  found  in  the  appendices  and  references[17-19]. 

4.1  ACCOMPLISHMENTS  OF  FIRST  TWO  YEARS 

The  first  two  years  of  Grant  AFOSR-91-0318  have  focused  on  selection  of  a 
numerical  scheme  that  (1)  is  appropriate  for  supercomputing,  (2)  that  scales  from  quasi- ID 
through  2-D  to  fully  3-D  problems,  and  (3)  which  can  be  easily  extended  to  hundreds  of 
states  associated  with  the  different  species  present  in  the  supersonic  nozzle  of  the  hydrazine 
and  ammonia  arcjets. 

Use  of  numerical  methods  for  highly  non-equilibrium  high  speed  reacting  internal 
flows  are  not  by  any  means  trivial.  Although  the  numerical  techniques  of  MacCormackf  20- 
22]  and  Beam  and  Warming[23]  are  well  developed  for  supersonic  external  compressible 
flows,  the  incorporation  of  finite  rate  kinetics  and  the  problem  of  internal  flow  complicates 
the  schemes  substantially.  Stringent  demands  must  therefore  be  placed  on  computational 
speed,  stability,  and  accuracy  especially  if  a  3-D  simulation  is  the  ultimate  goal.  The 
chosen  numerical  scheme  known  as  the  Linearized  Block  Implicit  (LBI)  method  first 
developed  by  Briley  and  McDonald[24]  shows  tremendous  promise.  It  is  the  only  method 
that  is  well  suited  to  supercomputers,  and  extends  all  the  way  from  quasi  1-D,  to  2-D  axi- 
symmetric,  through  to  fully  3-D  situations.  This  research  effort  has  been  divided 
into  2  two-year  increments  (overall,  a  four  year  effort),  in  which  the  first 
two  years  have  been  devoted  to  numerical  models  of  cold  flow,  ionizing 
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flow,  and  vibrationally  non-equilibrium  flow  under  quasi-ID,  2-D  with 
swirl,  and  fully  3-D  assumptions.  Each  of  these  cases  is  summarized  separately 
below,  but  extensive  details  may  be  found  in  the  appendices  as  well  as  in  the  annual  report 
for  the  first  year[25]. 

&  Quasi  1-P  Flow; 

Sample  calculations  of  quasi  1-D  flow  with  gas  internal-mode  molecular  energy  transfer 
in  the  30  KW  geometry  of  the  arc  jet  (see  Appendix  A)  have  been  performed  for  the  case 
of  smooth  acceleration  to  supersonic  speeds.  Each  of  these  energy  transfer  processes 
represent  elementary  mechanisms  and  could  just  as  easily  represent  chemical  reactions. 
Thus,  although  the  arc  and  propellant  chemistry  are  not  modelled  in  this  case,  this 
calculation  illustrates  the  power  of  the  numerical  algorithm  and  its  ultimate  usefulness 
for  simulating  real  arcjet  flows.  In  the  case  of  the  CO/N2  mixture  as  many  as  80  total 
internal  vibrational  states  are  modelled.  These  could  easily  have  been  as  many  chemical 
species,  or  represented  the  internal  vibrational  modes  in  a  N2/H2  mixture. 

(ii)  2-D  Axi-svmmetric  Flow  (with  swirll: 

Few  existing  works  have  successfully  examined  viscous  supersonic  internal  flows  with 
swirl.  We  have  completed  both  cold  flow  and  arcjet  simulations  (i.e.  governing 
equations  for  2-D  axi-symmetric,  viscous  flow  with  swirl,  combined  with  Maxwell’s 
equations,  and  with  finite  rate  ionization,  recombination,  and  dissociation  processes)  in 
the  30  KW  arcjet  geometry.  The  governing  equations  and  details  of  the  application  of 
the  LBI  method  to  these  equations  is  given  in  the  appendices.  Detailed  examination  of 
the  effect  of  swirl  on  cold  flow  in  geometries  of  interest  to  the  arcjet  thruster  reveal  that 
the  upstream  plenum  length,  i.e.  distance  from  the  swirl  injector  to  the  begining  of  the 
converging  section,  is  a  critical  design  parameter  for  a  given  operating 
propellant  mass  flow  rate.  For  the  case  of  high  mass  flow  rates  on  the  order  of 
g/s,  the  swirl  component  of  the  absolute  velocity  is  found  to  decrease  in  the  stream  wise 
direction  so  that  by  the  time  the  constrictor  is  reached  the  inlet  swirl  fails  to 
persist.  On  the  other  hand,  for  mass  flow  rates  on  the  order  of  mg/s,  the  swirl 
component  of  the  absolute  velocity  increases  in  the  streamwise  direction  and  reaches 
a  maximum  value  in  the  constrictor  region. 

The  enhancement  of  injected  swirl  may  be  explained  as  follows.  Injection  of  swirl 
into  the  arcjet  is  equivalent  to  an  influx  of  angular  momemtum  in  the  streamwise 
direction.  If  the  torque  due  to  the  shear  stress  at  the  wall  is  insufficient  to  overcome 
this  influx  of  angular  momentum,  then  the  inlet  swirl  will  persist  far  enough 
downstream  so  that  the  decreasing  cross-sectional  area  in  the  converging  section  of  the 
arcjet  geometry  will  cause  an  increase  in  the  swirl.  Although  this  result  is  observed  at 
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low  mass  flow  rates  (i.e.  on  the  order  of  mg/s),  it  strongly  suggests  that  for  any 
given  mass  flow  rate,  there  exists  a  maximum  plenum  length  beyond 
which  the  injected  swirl  will  have  no  effect  whatsover  on  the  arc 
downstream.  This  is  because  the  torque  due  to  the  wall  shear  stress  acting  on  the 
plenum  boundaries  is  large  enough  to  consume  the  inlet  flux  of  angular  momentum. 

Illustrative  results  for  arcjet  flow  are  shown  in  Figures  1  through  14.  These 
illustrative  cases  are  for  80  mg/s  of  hydrogen  flow  at  100  A  in  the  30  KW  geometry.  A 
full  discussion  of  these  and  other  results  for  argon  and  nitrogen  propellants  are 
discussed  in  a  forthcoming  paper[26]. 

(ui)  Fully  3-D  Simulations: 

Development  of  a  fully  3-D  cold-flow  code  (using  the  LBI  method)  from  a  2-D  axi- 
symmetric  cold  flow  model  is  relatively  straight  forward.  Although  such  a  3-D  code 
has  been  tested,  it  has  not  been  run  extensively  mainly  in  order  to  save  computer  time. 
We  envision  extensive  numerical  studies  being  conducted  using  our  2-D  axi-symmetric 
model,  and  the  fully  3-D  model  being  used  for  select  cases  once  the  intricate  details  of 
the  arcjet  flow  problem  have  been  fully  explored. 

The  fully  3-D  model  is  expected  to  be  useful  in  situations  where  conditions  of 
azimuthal  symmetry  break  down.  Such  a  situation  can  very  easily  occur  in  practice, 
where  the  centerline  of  the  cathode  and  anode  may  not  coincide  due  to  slight 
misalignment  during  initial  assembly.  The  resulting  effects  on  any  measurements  made 
downstream  of  the  cathode  tip  could  be  quite  severe.  This  will  be  further  explored. 

4.2  ON-GOING  &  FUTURE  WORK 

Research  has  been  initiated  in  four  parallel  areas,  namely  (1)  incorporation  of 
species  diffusion  (completed  as  of  March  1994),  (2)  incorporation  of  vibrational  and 
electronic  non-equilibrium  for  diatomic  species,  (3)  incorporation  of  governing 
equation  for  electron  energy  distribution  function,  and  (4)  modelling  of  anode 
attachment  without  ad-hoc  fixes.  Most  of  this  will  be  completed  in  the  second  two 
years  of  this  first  4-year  effort.  The  second  2-year  effort  will  involve  completion  of 
model  development,  extensive  comparison  with  available  experimental  data, 
suggestions  for  new  experiments,  and  detailed  performance  calculations  for  the  1  KW 
and  30  KW  arcjets. 

In  comparison  with  our  originally  proposed  work  statement,  we  are  exactly  on 
schedule  for  completion  of  the  proposed  tasks  for  the  second  two-years  of  the  full  4- 
year  period.  Since  our  approach  to  modelling  of  the  arcjet  flow  is  unique,  novel,  and 
extends  the  state-of-the-art,  we  anticipate  greater  interaction  with  experimentalists 
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making  measurements  on  the  1 KW  and  30  KW  axcjet  thrusters. 
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Fig.  I  :  Contours  of  y  (=  rB)  or  enclosed  current  contours  are  shown  here  for  the  hydrogen 

arcjet  at  100  A  and  a  mass  flow  rate  of  100  mg/s  in  the  30  KW  geometry.  Details  of  the 
calculation  are  given  in  the  appendices  as  well  as  in  a  forthcoming  papei{26]. 
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30  kW  Arcjet  Geometry 
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Fig.  2:  Contours  of  temperature  (non-dimensionalized  by  10,000  K)  are  shown  here  for  the  one- 
temperature  model.  The  inner-most  contour  represents  50,000  K,  while  the  outer-most 
represents  2,000  K. 
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Fig.  3:  Contours  of  the  Mach  number  (Flow  velocity  divided  by  the  frozen  speed  of  sound)  are 
shown  here.  M=1.3  is  the  right-most  contour 


Axial  Velocity 
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Fig.  4:  Contours  of  axial  component  of  velocity  (non-dimensionalized  by  8,323.6  m/s)  are 
shown  here.  Note  the  pockets  of  recirculating  flow  at  the  cathode  tip  as  represented  by 
negative  axial  velocity  contours.  These  pockets  have  recently  beer  observed 
experimentally  by  Cappelli  at  Stanford. 


versus  radius 


Axial  Velocity 


Fig.  6:  The  axial  velocity  (non-dimensionalized  by  8323.6  m/s)  is  shown  here  versus  radius 
(non-dimensionalized  by  0.01  m)  is  shown  here  at  the  exit  plane. 
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Fig.  8:  Profiles  of  temperature  (non-dimensionalized  by  10,000  K)  are  shown  here  versus  x 
(non-dimensionalized  by  0.01  m)  at  various  radial  locations.  j=0  corresponds  to  the 
cathode  surface  and  channel  centerline,  while  j=60  corresponds  to  locations  near  the 
anode. 
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Fig.  9: 


Profiles  of  axial  component  of  velocity  (non-dimensionalized  by  8323.6  m/s)  are 

^r^  ^anode.  Note  dte  presence  of  d«  mcircnladon  zones  or  separaoon 

bubbles  at  the  cathode  tip. 


Temperature 


Fig.  10:  Profiles  of  temperature  (non-dimensionalized  by  10,000  K)  are  shown  here  versus  r 
(non-dimensionalized  by  0.01  m)  at  various  axial  locations.  i=45  corresponds  to 
locations  near  the  cathode  dp,  while  i=70  corresponds  to  locations  after  the  constrictor 
and  in  the  nozzle.  i=50,  i=55,  and  i=60  all  represent  locations  near  the  constrictor  inlet, 
constrictor  center,  and  constrictor  exit  respectively. 
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Axial  Velocity 


Fig.  11:  Profiles  of  axial  component  of  velocity  (non-dimensionalized  by  8323.6  m/s)  are 
shown  here  versus  r  (non-dimensionalized  by  0.01  m)  at  various  axial  locations.  i=45 
corresponds  to  locations  near  the  cathode  tip,  while  i=70  corresponds  to  locations  after 
the  constrictor  and  in  the  nozzle.  i=50,  i=55,  and  i=60  all  represent  locations  near  the 
constrictor  inlet,  constrictor  center,  and  constrictor  exit  respectively. 


Electron  Number 

Fig.  12:  Profiles  of  electron  number  density  (non-dimensionalized  by  5.8X1022  m'3)  are  shown 
here  versus  r  (non-dimensionalized  by  0.01  m)  at  various  axial  locations.  i=45 
corresponds  to  locations  near  the  cathode  tip,  while  i=70  corresponds  to  locations  after 
the  constrictor  and  in  the  nozzle.  i=50,  i=55,  and  i=60  all  represent  locations  near  the 
constrictor  inlet,  constrictor  center,  and  constrictor  exit  respectively.  Note  the  sharp 
gradients  in  the  profiles  which  clearly  define  the  arc  boundary  in  the  absence  of 
diffusion.  The  electron  number  density  is  seen  to  increase  with  radius  and  then 
abruptly  decrease  outside  the  arc  boundary.  Such  maxima  are  typical  when  there  are 
opposing  mechanisms,  i.e  ohmic  heating  tending  to  increase  ionization  while  the 
temperature  decreasing  with  radius  tends  to  decrease  the  ionization. 


H-atom  number 


Fig.  13:  Profiles  of  atomic  hydrogen  number  density  (non-dimensionalized  by  5.8X1022  m"3) 
are  shown  here  versus  r  (non-dimensionalized  by  0.01  m)  at  various  axial  locations. 
i=45  corresponds  to  locations  near  the  cathode  tip,  while  i— 70  corresponds  to  locations 
after  die  constrictor  and  in  the  nozzle.  i=50,  i=55,  and  i=60  all  represent  locations  near 
the  constrictor  inlet,  constrictor  center,  and  constrictor  exit  respectively.  Note  that  here 
in  the  absence  of  diffusion,  the  atomic  hydrogen  concentration  has  a  maxima  in  the 
same  location  as  the  electron  number  density. 
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Fig.  14:  Profiles  of  molecular  hydrogen  number  density  (non-dimensionalized  by  5.8x1 022  m  ) 
are  shown  here  versus  r  (non-dimensionalized  by  0.01  m)  at  various  axial  locations. 
i=45  corresponds  to  locations  near  the  cathode  tip,  while  i=70  corresponds  to  locations 
after  the  constrictor  and  in  the  nozzle.  i=50,  i=55,  and  i=60  all  represent  locations  near 
the  constrictor  inlet,  constrictor  center,  and  constrictor  exit  respectively.  Note  that  the 
molecular  hydrogen  concentration  is  largest  near  the  anode,  where  both  electron  and 
atomic  hydrogen  concentrations  are  smaller. 
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1  Introduction 

The  compressible  Navier-Stokes  equations  govern  the  physics  of  flow  fields  in  many  applications. 
Recently,  many  aerospace  applications  involving  both  internal  and  external  supersonic  flows 
have  arisen.  These  flows  are  often  times  further  complicated  by  the  effects  of  departure  from 
thermodynamic  equilibrium  of  the  internal  modes  of  molecular  motion.  This  disequilibrium 
between  internal  modes  (vibration,  rotation  and  electronic  excitation)  and  external  modes  i.e. 
translational  molecular  motions  gives  rise  to  complicated  chemistry.  Before  such  complex  flows 
can  be  studied  however,  the  basic  supersonic  viscous  flow  equations  must  be  solved  using  a 
reliable  numerical  method  that  can  be  readily  extended  to  multi-dimensions.  This  is  the  focus 
of  the  present  article. 

Many  numerical  methods  have  been  devised  to  solve  the  compressible,  Navier-Stokes  equa¬ 
tions,  the  most  popular  of  these  being  the  MacCormack  scheme  [1],  the  Beam- Warming  scheme 
[2]  and  the  linearized  block  implicit  (LBI)  scheme  of  Briley  &  McDonald  [3].  There  are  many 
solutions  for  viscous,  external,  supersonic  flows  obtained  with  any  of  these  methods  (see  for 
example  [4,  5]),  and  there  are  a  few  recent  works  that  include  simple  chemical  kinetics  [6,  7).  In 
contrast,  there  are  only  a  limited  number  of  available  numerical  solutions  for  internal,  viscous 
supersonic  flows  [8].  Further,  there  has  been  only  one  recent  work  that  investigates  the  effects  of 
non-equilibrium  state  resolved  chemistry  and  vibrational  non-equilibrium  in  a  supersonic  expan¬ 
sion  in  a  converging-diverging  (C-D)  nozzle  [9].  This  work  however  examines  the  flow  under  the 
quasi  one-dimensional  assumption.  Our  ultimate  goal  is  to  obtain  two-  and  three-dimensional 
solutions  that  include  such  state  resolved  kinetics. 

In  this  paper,  numerical  solutions  to  viscous  subsonic  and  supersonic  flows  in  C-D  nozzles 
are  presented  both  with  and  without  the  effects  of  swirl.  The  area  ratios  of  the  nozzles  that  we 
consider  are  high,  being  of  the  order  of  25  and  can  be  even  higher.  The  LBI  scheme  of  Briley 
&  McDonald  [3]  is  used  to  obtain  these  solutions.  The  importance  of  studying  such  flows  in 
this  geometry  cannot  be  emphasized  enough,  for,  the  essential  physics  of  internal,  supersonic 
flows  in  many  problems  of  practical  interest  can  be  related  to  the  flow  in  a  corresponding  C-D 
nozzle  geometry.  With  this  in  mind,  we  present  both  inviscid  quasi  one-dimensional  solutions 
and  two-dimensional  axisymmetric,  viscous  flow  solutions.  We  find  that  the  key  to  obtaining 
these  solutions  is  the  proper  modelling  of  the  exit  boundary  conditions  and  specification  of  the 
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proper  initial  conditions  so  that  the  flow  can  accelerate  smoothly  from  subsonic  to  supersonic 
speeds.  The  algorithm  itself  is  easily  extendable  to  three  dimensions. 

This  paper  is  organized  as  follows.  The  governing  equations  and  boundary  conditions  are 
presented  first,  followed  by  the  numerical  method.  Quasi  1-D  and  2-D  axisymmetric  solutions 
with  and  without  swirl  are  discussed  next,  followed  by  a  summary  and  conclusions  of  this  work. 


2  Formulation 


The  equations  that  govern  the  cold  (i.e.  non-reacting)  viscous,  supersonic  flow  axe  the  com¬ 
pressible  Navier-Stokes  equations.  For  our  present  purposes,  we  consider  them  in  cylindrical 
polar  coordinates  in  2-D  axisymmetric  form.  These  equations  are  made  dimensionless  by  non- 
dimensionalizing  coordinates  by  the  throat  radius,  density  by  the  inlet  stagnation  density,  tem¬ 
perature  by  the  inlet  stagnation  temperature,  and  velocities  by  the  inlet  stagnation  speed  of 
sound  respectively.  Pressure  is  eliminated  from  the  system  of  equations  by  using  the  equation  of 
state.  This  choice  of  non-dimensionalization  differs  from  the  conventional  ones  (5,  section  5-1.6] 
but  is  closely  related  to  the  one  used  in  [3].  Since  in  general  in  a  reacting  flow  the  properties  such 
as  specific  heats  and  viscosity  are  locally  varying  (and  in  most  instances  by  as  much  as  orders 
of  magnitude),  non-dimensionalization  using  quantities  such  as  local  frozen  speed  of  sound  [10], 
pressure  etc.,  are  not  useful.  We  choose  therefore  as  reference  quantities,  stagnation  conditions 
upstream  at  the  inlet  and  as  will  be  seen  later,  two  of  these  quantities,  namely,  stagnation 
pressure  and  stagnation  temperature  have  to  be  prescribed  at  the  inlet  as  boundary  conditions. 
Accordingly,  we  have  in  conservative,  non-dimensional  form: 
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where,  $  is  the  dissipation  function  given  by, 
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where,  A  —  Arcr  +  Age's  +  Aze~x  is  any  vector  and  <f>  is  a  scalar. 


In  the  above  equations,  tT  =  ue*r  +  vt»  +  toe',  is  the  velocity  vector,  7  =  Cp/Cv  is  the 
ratio  of  specific  heats,  Re  =  pa<ioUhToa.tl  P  is  the  Reynolds  number,  and  Pr  =  pCv/k  is  the 
Prandtl  number.  The  subscript  o  denotes  stagnation  conditions  at  the  inlet  and  a0  =  y/f  RT0 
is  the  stagnation  speed  of  sound  at  the  inlet.  It  is  assumed  throughout  in  this  paper  that  the 
coefficient  of  viscosity  p,  the  thermal  conductivity  k  of  the  fluid,  together  with  other  properties 
such  as  7  remain  constant  at  representative  values.  Variable  properties  can  be  easily  handled 
by  the  LBI  scheme  [3]. 


2.1  Boundary  Conditions 

Boundary  conditions  for  the  five  dependent  variables  (p,u,v,w,  T)  need  to  be  specified  along 
the  inlet  and  exit  boundaries,  on  the  wall,  and  along  the  centerline  of  the  nozzle  geometry. 
No-slip  conditions  are  imposed  on  the  wall  which  is  maintained  at  a  constant  temperature.  The 
density  at  the  walls  can  be  specified  either  by  applying  the  unsteady  continuity  equation  [3]  or 
the  normal  component  of  the  momentum  equation  on  the  wall.  We  choose  the  latter,  since  it 
is  physically  more  meaningful  to  cast  the  boundary  conditions  in  terms  of  momentum  flux,  but 
apply  the  radial  component  of  the  momentum  equation  instead  of  the  normal  component  on  the 
wall.  The  choice  of  the  radial  component  instead  of  the  normal  component  is  dictated  by  the 
need  for  simplicity.  Thus  we  have. 


u  a  V  =  W  =  0,  T  =  T,urfaee,  ©  r  =  R(z), 
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where  r  =  R(z)  is  the  equation  for  the  outer  wall  of  the  nozzle. 

Symmetry  conditions  are  imposed  along  the  centerline  for  all  the  variables  except  the  swirl 
component  of  velocity  t>,  which  is  set  equal  to  zero,  identically.  These  are  written  as, 
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Prescribing  boundary  conditions  at  the  inlet  and  exit  planes  of  the  nozzle  depends  on  whether 
the  flow  is  subsonic  or  supersonic.  All  the  solutions  presented  here  consider  subsonic  inlets. 
Hence,  following  Briley  &  McDonald  [3],  the  stagnation  temperature  and  stagnation  pressure 
at  the  inlet  are  specified.  In  addition,  at  the  inlet,  the  radial  component  of  velocity  is  set  to 
zero  and  the  swirl  component  of  velocity  is  either  zero  or  set  to  some  prescribed  distribution 
corresponding  to  injection.  It  is  important  to  note  here  that  the  axial  component  of  velocity 
which  determines  the  mass  flow  rate  through  the  nozzle,  is  not  specified  at  the  inlet  but  is 
determined  from  an  extrapolation  of  the  interior  values  (see  the  difference  equation  (11)  below). 
The  above  conditions  can  be  written  in  dimensionless  form  as, 

u  =  0,  v  =  V(r),  @  z  =  0,  (9) 

T+  +  1,  =  1,  @z  =  0.  (10) 

^  y  t— r 

Along  the  exit  plane,  the  flow  is  entirely  subsonic  in  the  case  of  converging  nozzles  or  mixed 
subsonic-supersonic  in  the  case  of  the  C-D  nozzles.  In  either  case,  all  the  variables  except  density 
on  the  exit  plane  axe  evaluated  by  extrapolation  from  the  interior  [3].  If  the  local  Mach  number 
AT  defined  as  t o/y/T  is  greater  than  one  at  a  point  on  the  exit  plane,  then  the  density  is  also 

evaluated  by  extrapolation.  If,  on  the  other  hand,  M  <  1  at  some  point  on  the  exit  plane,  then 

the  exit  static  pressure  is  specified  and  the  density  is  evaluated  from  the  equation  of  state.  In 
the  case  of  C-D  nozzles,  the  subsonic  flow  occurs  within  the  boundary  layer  near  the  wall,  and 
the  correct  value  of  the  exit  static  pressure  would  usually  not  be  known,  unless  the  plume  is 
modelled  as  well.  Care  must  be  exercised  in  choosing  a  value  for  the  exit  static  pressure,  for,  if 
it  is  too  high,  the  boundary  layer  on  the  wall  in  the  diverging  section  of  the  nozzle  can  separate. 
This  is  an  important  detail  that  often  escapes  mention  in  the  literature. 

3  Numerical  Method 

As  a  first  step  towards  solving  equations  (l)-(5),  a  non-staggered  (non-uniform  or  uniform)  grid 
is  generated  in  the  physical  domain.  The  non-uniform  grids  are  generated  using  the  Roberts 
transformation  [5,  section  5-6.1].  This  grid  in  the  physical  domain  is  then  transformed  into  a 
uniform  grid  in  a  computational  domain  using  an  analytical  coordinate  transformation.  The 
transformation  is  such  that  the  nozzle  contour  is  transformed  into  a  rectangle  in  the  compu¬ 
tational  domain  [5,  section  5-6.1].  This  transformation  can  be  conveniently  written  as  follows 
[12]: 

r=f  +  |/?tanh[(^-l/2)ln(|±|)]  0  <  r  <  R,  0  <  f  <  1, 

(  =  z, 

where,  {  is  the  radial  coordinate  and  (  is  the  axial  coordinate  in  the  computational  domain. 
The  quantity  0  is  a  stretching  parameter.  It  is  set  equal  to  1.5  for  generating  the  non-uniform 
grid  mentioned  later  in  Section  4.  In  the  limit  as  0  -*  oo  the  above  transformation  yields  a 
uniform  grid  in  the  radial  direction  also.  Equations  (l)-(5)  together  with  the  boundary  conditions 
(6)-(ll)  are  then  transformed  to  the  computational  domain  and  solved. 
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The  transformed  forms  of  these  equations  are  solved  using  the  LB1  scheme  of  Briley  & 
McDonald  (11).  Since  all  the  details  of  the  scheme  have  already  been  reported  elsewhere  [11,  3], 
only  the  salient  features  will  be  discussed  here.  First,  the  equations  are  discretized  in  time  using 
Crank- Nicolson  differencing.  Even  though  this  scheme  is  second  order  accurate,  the  linearization 
of  the  non-linear  time  derivatives  is  only  first  order  accurate,  so  that  the  algorithm  overall  is 
only  first  order  accurate  in  time.  This  is  not  a  major  concern  if  only  steady  state  solutions  are 
desired,  as  is  the  case  here.  The  non-linear  spatial  derivatives,  however,  are  linearized  to  second 
order  accuracy  in  time.  These  and  the  other  spatial  derivatives  of  the  dependent  variables  are 
then  discretized  using  second  order  accurate  central  differences.  The  transformation  derivatives 
are  evaluated  analytically.  The  boundary  conditions  also  involve  non-linear  quantities  and  are 
linearized  to  second  order  time  accuracy  in  the  same  manner  as  the  governing  equations.  All 
the  spatial  derivatives  that  appear  in  these  boundary  conditions  are  evaluated  using  second 
order  accurate  one-sided  differences.  The  extrapolation  that  is  required  in  the  inlet  and  exit 
boundaries  can  be  done  implicitly  in  a  manner  suggested  by  Briley  &  McDonald  [3].  This  can 
be  written  in  difference  form  at  the  inlet,  for  example,  as, 

V#1  =  2V#1  -  V#1  1  <  3  <  M  ,  (11) 

where,  ip  is  any  variable  that  needs  to  be  extrapolated.  The  superscript  n+1  denotes  the  new 
time  level  and  the  subscripts  denote  the  grid  indices  in  the  usual  sense.  Similar  expressions  are 
written  for  the  exit  plane. 

The  resulting  coupled,  linear  algebraic  equations  are  then  solved  by  using  the  Douglas-Gunn 
ADI  method.  If  this  system  of  equations  were  written  in  matrix  form,  then  solution  by  the  ADI 
method  can  be  viewed  as  computing  the  inverse  of  the  matrix  to  within  a  certain  splitting  error 
which  is  of  the  same  order  (or  higher)  as  the  time  discretization  error  in  the  original  equations 
[11].  The  ADI  method  itself  involves  the  solution  of  block  tridiagonal  systems  in  each  coordinate 
direction,  with  block  size  equal  to  the  number  of  dependent  variables,  which,  in  our  case  is  five. 
These  systems  can  be  solved  using  existing  efficient  LU  decomposition  methods  [3]. 


3.1  Artificial  Dissipation 


It  is  well  known  that  artificial  dissipation  of  some  form  is  necessary  [4,  section  18.5]  in  order 
to  avoid  spurious  oscillations  (and  possible  loss  of  convergence)  especially  at  high  Reynolds 
numbers.  This  difficulty  arises  due  to  the  central  differencing  [13,  section  III-C-8]  used  for 
discretizing  the  convective  terms  in  equations  (l)-(5).  Roache  [13]  shows  that  numerical  solutions 
exhibit  non-physical  behavior  if  the  cell  Reynolds  number  Re&x  =J  to  |  Az  Re  is  greater  than 
two.  This  numerical  instability  can  be  avoided  by  using  upwinded  differencing  for  the  convective 
terms  or  by  adding  explicitly  a  certain  amount  of  second  order  or  fourth  order  dissipation  [4, 
p.439].  Briley  &  McDonald  [3]  used  a  second  order  dissipation  term  of  the  form  where 


e*  = 


to 


Az 


-  >  2, 


€*  =  0,  Re&t  <  2. 

This  assures  that  the  local  cell  Reynolds  number  is  less  than  or  equal  to  two.  This  dissipation 
term  is  added  to  each  one  of  the  governing  equations  with  <f>  defined  approriately  as  in  [3].  This 
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is  done  in  the  present  work  also.  The  dissipation  term,  while  improving  the  stability  of  the 
numerical  method,  has  a  detrimental  effect  on  the  accuracy  of  the  solution  [3].  It  is  difficult 
to  ascertain  this  loss  of  accuracy  in  a  multidimensional  problem  such  as  the  one  we  are  dealing 
with  here.  Nevertheless,  an  attempt  was  made  to  minimize  this  loss  of  accuracy  through  the 
use  of  an  adequate  number  of  grid  points  such  that  all  the  scales  of  interest  in  the  flow  field  are 
resolved  properly. 

4  Results  &  Discussion 

In  this  section,  we  first  present  quasi  1-D  solutions  for  the  flow  in  a  C-D  nozzle,  followed  by 
viscous  2-D  axi-symmetric  solutions.  The  quasi  I-D  solutions  are  obtained  for  the  flow  of  air 
through  a  nozzle  whose  area  variation  is  given  as  A(z)  =  1  +  z2,  —  1  <  z  <  10.  This  profile 
gives  an  exit  to  throat  area  ratio  of  approximately  100.  Two  solutions  are  presented  in  Figure 
1  -  one  with  an  exit  pressure  high  enough  to  initiate  a  normal  shock  in  the  diverging  section 
of  the  nozzle,  and  another  where  the  flow  accelerates  smoothly  from  a  subsonic  velocity  at  the 
inlet  to  supersonic  velocities  at  the  exit.  The  fully  supersonic  solution  has  been  obtained  with 
zero  artificial  dissipation,  while  the  other  solution  has  been  obtained  with  an  amount  of  artificial 
dissipation  sufficient  to  make  the  local  cell  Reynolds  number  equal  to  2.  A  linearly  decreasing 
density  profile  and  constant  values  for  the  velocity  (typically  zero)  and  temperature  (typically 
278  K)  are  used  as  initial  guesses  to  start  the  solution  procedure.  A  total  of  256  points  are  used 
and  the  calculations  are  stopped  when  the  percent  change  in  the  dependent  variables  is  less  than 
0.00001%.  As  a  check  on  the  accuracy  of  the  solution  the  mass  flow  rate  is  computed  at  each 
axial  station  in  the  nozzle,  and  is  found  to  be  constant  to  within  0.00001%.  These  solutions 
take  about  10  to  15  seconds  of  CPU  time  on  the  Cray  Y-MP.  It  can  be  seen  in  Figure  1  that 
the  two  curves  are  barely  distinguishable  from  each  other  up  until  the  shock,  which  shows  that 
the  effect  of  the  added  artificial  dissipation  on  the  accuracy  of  the  solution  is  minimal.  We  now 
turn  to  a  discussion  of  the  viscous  flow  solutions. 

All  the  viscous  flow  solutions  presented  here  are  for  the  flow  of  air  (7  =  1.4,  Pr  =  0.701) 
through  a  nozzle  whose  geometry  is  shown  in  Figure  2.  It  has  a  converging  and  a  diverging 
section  in  which  the  outer  radius  of  the  nozzle  is  a  linear  function  of  the  axial  coordinate.  The 
transition  from  the  converging  to  the  diverging  section  is  achieved  within  five  grid  points  through 
the  use  of  a  quadratically  varying  profile  for  the  outer  radius.  The  exit  to  throat  area  ratio  for 
this  nozzle  geometry  is  25.  Our  interest  in  this  geometry  arises  from  its  use  in  space  propulsion 
applications  [14].  All  the  solutions  presented  here  are  obtained  on  a  non-uniform  grid  with  256 
points  in  the  axial  direction  and  64  points  in  the  radial  direction.  The  calculations  are  performed 
until  the  percent  change  in  the  dependent  variables  is  less  than  0.00001%. 

It  is  essential  to  specify  proper  initial  conditions  in  order  to  obtain  the  supersonic  solutions 
that  are  presented  here.  While  this  may  not  be  necessary  in  the  case  of  C-D  nozzles  with  smaller 
area  ratios,  we  find  this  to  be  crucial  for  the  large  area  ratio  that  we  consider  here.  The  key 
feature  of  this  initial  guess  is  that  it  must  allow  the  C-D  nozzle  to  choke  very  quickly.  This  is 
important  because  once  the  nozzle  chokes,  the  flow  in  the  diverging  section  reaches  supersonic 
velocities  quite  rapidly.  The  appropriate  initial  guess  can  be  generated  by  solving  the  flow  in  a 
geometry  that  has  the  same  converging  section  as  the  C-D  nozzle  followed  by  a  constant  area 
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section  (with  area  equal  to  the  throat  area)  in  place  of  the  diverging  section.  The  exit  static 
pressure  of  this  geometry  is  adjusted  until  a  converged  solution  with  a  maximum  exit  Mach 
number  close  to  but  less  than  one  is  obtained.  This  solution  is  then  used  as  the  initial  guess  for 
generating  the  flowfield  in  the  C-D  nozzle.  For  the  present  case,  on  the  256x64  grid,  it  takes 
approximately  4-5  minutes  of  CPU  time  on  a  Cray  Y-MP  to  generate  this  initial  guess.  The 
desirable  feature  of  this  procedure  is  that  it  is  easy  to  generate  this  initial  guess  even  in  complex 
geometries  with  minimal  computational  effort.  With  this  initial  guess,  it  takes  approximately 
8-10  minutes  of  CPU  time  on  the  Cray  to  obtain  the  solutions  for  the  C-D  nozzle,  which  are 
discussed  in  the  following  paragraphs. 

The  velocity  vectors  for  the  flow  in  the  C-D  nozzle  at  a  Reynolds  number  (as  defined  in 
Section  2)  of  50,000  are  shown  in  Figure  3.  A  throat  radius  of  2  mm,  inlet  stagnation  pressure 
of  1.03  x  10s  N/m2,  inlet  stagnation  temperature  of  300  K  and  a  coefficient  of  viscosity  p  = 
18  x  10~6  Ns/m 2  are  used  for  calculating  this  Reynolds  number.  Different  values  of  Reynolds 
numbers  can  be  obtained,  for  example,  by  adjusting  the  inlet  stagnation  quantities.  The  velocity 
vectors  show  quite  clearly  the  growth  of  the  boundary  layer  near  the  wall  in  the  diverging  section. 
Only  one  half  of  the  total  number  of  points  in  the  radial  direction  and  one  tenth  of  the  total 
number  in  axial  direction  are  shown  in  this  figure  for  the  sake  of  clarity.  One  can  also  see  some 
wall- jet  type  velocity  profiles  in  the  vicinity  of  the  throat  where  the  flow  accelerates  rapidly  in 
order  to  negotiate  the  sharp  corner. 

Profiles  of  Mach  number  along  the  radial  grid  lines  j  =  0  and  j  =  38  are  shown  in  Figure  4 
together  with  the  Mach  number  profile  from  the  corresponding  frictionless,  quasi  1-D  solution. 
Here,  the  grid  line  corresponding  to  j  =  0  represents  the  centerline.  The  profile  along  the  j  =  38 
grid  line  is  closest  to  the  quasi  1-D  profile.  The  maximum  exit  Mach  number  in  the  viscous  case 
is  about  4.6,  which  is  less  than  the  value  of  about  5.0  for  the  frictionless,  quasi  1-D  case.  The 
profile  along  the  centerline  shows  a  sharp  drop  in  the  Mach  number  near  a  dimensionless  axial 
coordinate  of  8.  One  can  see  from  the  Mach  number  contours  in  Figure  5,  a  change  in  the  shapes 
of  these  contours  near  the  same  axial  location  (z  =  8).  In  this  figure,  the  sonic  line  (Af=l)  is 
shown  as  the  thick  line.  The  radial  Mach  number  profiles  have  their  maxima  located  on  the 
centerline  up  until  z  =  8  and  these  maxima  move  to  a  location  near  the  center  of  the  channel  once 
we  cross  this  particular  axial  location.  These  contours  look  very  similar  to  those  obtained  using 
the  method  of  characteristics  [15]  for  inviscid  flow  in  a  similar  geometry.  They  also  agree  very 
well  qualitatively  with  the  experimental  results  presented  in  [16].  It  was  first  shown  by  Darwell 
&  Badham  [15]  that  compression  waves  could  be  generated  in  C-D  nozzle  geometries  from  the 
location  where  the  curved  throat  section  meets  the  diverging  conical  section,  if  the  radius  of 
curvature  of  the  nozzle  wall  profile  is  discontinuous  at  this  location.  These  compression  waves 
can  coalesce  at  the  centerline  at  a  point  downstream  to  form  a  weak  oblique  shock.  Experiments 
conducted  subsequently  by  Back  et  al.  [17,  18]  and  Cuffel  et  al.  [16]  showed  this  to  be  indeed 
true.  The  sharp  drop  in  the  Mach  number  along  the  centerline  in  Figure  4  is  due  the  oblique 
shocks  intersecting  the  centerline  and  the  change  in  the  shapes  of  the  Mach  number  contours 
mentioned  before  is  due  to  these  incoming  and  reflected  weak  oblique  shock  waves.  Indeed,  the 
similarity  between  Figure  4  and  the  experimental  profiles  of  Back  et  al.  [17,  Figure  2]  is  striking. 

The  exact  mechanism  by  which  these  compression  waves  are  generated  is  discussed  in  detail 
in  [15,  18].  As  the  flow  passes  through  the  throat,  the  regions  of  the  flow  close  to  the  wall 
are  overturned  so  as  to  be  parallel  to  the  wall.  Consequently,  when  the  location  where  the 


7 


discontinuous  change  in  wall  curvature  occurs  is  reached,  the  flow  must  be  turned  back  towards 
the  centerline,  i.e.  away  from  the  wall.  This  can  only  occur  through  a  series  of  compression 
waves.  This  can  be  seen  quite  clearly  in  Figure  6,  where  the  dimensionless  wall  static  pressure 
is  plotted.  A  bump  can  be  seen  in  the  curve  immediately  past  the  the  throat  section  which  is 
located  at  z= 4,  This  bump  is  shown  enlarged  in  the  inset.  Again,  this  profile  shows  excellent 
qualitative  agreement  with  the  corresponding  plot  in  [18,  Figures  3,5]. 

The  existence  of  the  incident  and  the  reflected  oblique  shock  waves  can  be  seen  clearly  in 
the  next  two  figures.  Figure  7  shows  a  plot  of  the  radial  velocity  along  the  grid  line  j  =  8. 
The  radial  velocity  is  negative  in  the  converging  section  and  then  it  rapidly  turns  around  at  the 
throat  becoming  positive  in  the  diverging  section.  The  radial  velocity  then  decreases  sharply 
in  the  neighbourhood  of  z  =  8  and  becomes  negative  as  the  flow  traverses  the  incident  oblique 
shock  waves.  The  fluid  is  turned  in  a  clockwise  direction  towards  the  centerline  by  these  waves, 
since  they  are  right  running.  The  radial  velocity  then  increases  just  as  rapidly,  as  the  flow 
traverses  the  oblique  shocks  that  are  reflected  from  the  centerline.  These  turn  the  fluid  in  an 
anti-clockwise  direction  away  from  the  centerline,  being  left  running  waves.  The  radial  velocity 
settles  down  to  a  constant  value  after  about  z  =  14.  The  overturning  of  the  fluid  as  it  goes 
through  the  throat  is  evident  from  this  figure  when  one  compares  the  value  of  the  radial  velocity 
near  z  «  7  with  the  final  value  for  z  >  14. 

The  total  pressure  along  the  centerline  (j  =  0)  and  along  the  radial  grid  line  j  =  10  is 
plotted  in  Figure  8.  The  total  pressure  is  essentially  constant  in  the  converging  section  where 
the  boundary  layer  is  developing  so  that  the  flow  is  mostly  inviscid.  But  viscous  effects  become 
important  in  the  throat  region  and  the  total  pressure  increases  as  a  result  of  heating  due  to 
viscous  dissipation.  The  total  pressure  then  decreases  continuously  as  the  flow  accelerates  to 
supersonic  speeds  and  increases  across  the  incident  and  reflected  oblique  shock  waves.  In  the 
case  of  the  profile  along  the  centerline,  this  increase  is  quite  a  sharp  one,  but  away  from  the 
centerline,  along  j  =  10,  two  bumps  in  the  profile  can  be  seen,  each  representing  the  pressure  rise 
following  the  incident  and  reflected  oblique  shocks  respectively.  This  agrees  quite  well  with  the 
experimental  results  ([17,  Figure  1].  Mach  number  contours  (Figure  9)  and  profiles  (Figure  10) 
for  a  Reynolds  number  of  1000  (corresponding  to  a  lower  mass  flow  rate)  show  that  the  effects 
of  the  rapid  turning  near  the  throat  are  mitigated  in  a  fully  viscous  flow.  This  is  to  be  expected 
since  at  lower  mass  flow  rates  and  lower  inlet  stagnation  pressures,  the  viscous  effects  dominate 
the  flow. 

We  now  consider  briefly  the  case  where  the  swirl  velocity  at  the  inlet  is  as  shown  in  Figure 
11.  This  type  of  swirl  profile  is  observed,  for  example,  if  the  fluid  were  injected  through  ports 
peripheral  to  the  nozzle  wall.  Thus  the  radial  velocity  at  inlet  is  still  zero,  as  before.  The 
flow  field  in  this  case  is  very  similar  to  the  previous  case  with  zero  swirl  and  so  they  are  not 
presented  here.  Contours  of  swirl  velocity  are  shown  in  Figure  12.  The  injected  swirl  can  be 
seen  to  decay  within  a  short  distance  from  the  entrance  due  to  viscous  effects.  However,  we  have 
observed  that  under  some  circumstances  this  injected  swirl  can  persist  and  can  get  enhanced 
quite  dramatically  in  geometries  of  the  type  used  in  [14]. 
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5  Summary  &  Conclusions 

Solutions  to  the  2-D,  axi-symmetric,  viscous,  compressible  flow  in  a  C-D  nozzle  with  a  large  exit 
to  throat  area  ratio  are  presented.  In  addition,  quasi  1-D  solution  for  flow  in  a  C-D  nozzle  with 
an  exit  to  throat  area  ratio  equal  to  100  is  also  given.  These  solutions  are  obtained  by  using  the 
LBI  scheme  of  Briley  &  McDonald  [3].  In  order  to  successfully  obtain  solutions  to  this  problem, 
proper  initial  and  boundary  conditions  need  to  be  specified  as  discussed  in  detail  in  this  paper. 
We  obtained  solutions  to  the  2-D,  axi-symmetric,  in  viscid  Euler  equations,  though  they  are  not 
presented  here.  This  system  of  equations  requires  artificial  dissipation  in  both  the  radial  and 
the  axial  directions  for  numerical  stability.  For  this  reason,  the  so  called  "inviscid”  solution  is 
not  truly  inviscid  but  resembles  the  viscous  solution  quite  closely.  Furthermore,  the  computing 
effort  required  is  approximately  the  same  in  both  the  cases  and  so  it  is  better  to  compute  the 
viscous  solutions  which  are  more  reflective  of  reality  than  the  corresponding  inviscid  solutions. 
Thus  the  quasi  1-D  solution  is  the  only  true  inviscid  solution. 

The  viscous  solution  shows  considerable  departure  from  the  quasi  1-D  solution.  Most  impor¬ 
tantly,  oblique  shocks  generated  due  to  overturning  of  the  fluid  dose  to  the  wall  in  the  throat 
region  are  seen  in  the  viscous  solution.  These  oblique  shocks  were  first  predicted  by  Darwell 

6  Badham  [15]  based  on  invisdd  flow  calculations  and  subsequently  confirmed  by  the  experi¬ 
ments  of  Back  et  al.  [17,  18].  The  qualitative  agreement  between  our  numerical  results  and  the 
experimental  results  are  excellent,  thus  illustrating  the  capability  of  the  LBI  scheme  to  capture 
the  true  physics  of  the  problem.  Addition  of  swirl  to  this  flow  does  not  produce  any  significant 
change  in  the  flow  field  and  the  injected  swirl  is  seen  to  decay. 

The  LBI  scheme  can  be  easily  extended  to  three  dimensions  and  can  account  for  variable 
properties  as  well.  It  retains  its  unconditional  stability  properties  in  three  dimensions  also  [11]. 
The  scheme  can  also  be  used  to  study  chemically  reacting  flows  and  its  capability  to  handle  the 
stiff  systems  that  arise  as  a  result  of  the  chemical  reactions  has  already  been  established  [7]. 
Further,  reacting  flows  such  as  those  studied  in  [9]  with  state  resolved  kinetics  can  be  effectively 
analyzed  using  the  LBI  method. 
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ABSTRACT 

Arcjets  operating  on  molecular  propellants 
such  as  ammonia  or  hydrazine  are  expected  to 
produce  significant  amounts  of  diatomic  species 
such  as  N2,  NH,  and  Hi  in  the  discharge 
[1].  Internal  mode  energy  exchange  (Vibration- 
Vibration,  Vibration- Electronic,  Electronic- 
Electronic,  and  Vibration-Translation)  be¬ 
tween  these  molecules  can  result  in  highly  non- 
Maxwellian  distribution  of  energies  within  their 
internal  modes  of  vibration,  electronic  excita¬ 
tion,  and  rotation.  These  diatomic  species  are 
also  extremely  effective  in  moderating  electron 
energies  via  inelastic  collisions  involving  their 
internal  modes  (vibration,  electronic  excita¬ 
tion,  and  rotation).  Consequently  the  electron 
energies  in  arcjets  operating  on  such  molecu¬ 
lar  propellants  are  strongly  related  to  and  de¬ 
termined  by  the  distribution  of  molecules  in 
various  vib»-'ionally  and  electronically  excited 
states.  U  r  these  conditions  therefore,  the 
concept  o*  an  electron  temperature  loses  mean¬ 
ing.  Instead,  one  must  refer  to  electron  energy 
distribution,  or  an  average  electron  energy  if 
a  meaningful  average  exists.  The  modelling  of 
such  inter-mode  energy  transfer  is  extremely 
complicated  and  presents  a  challenge  to  the  ar¬ 
eas  of  chemical  physics  and  numerical  analysis. 
In  this  paper,  we  present  numerical  solutions 
to  quasi  one-dimensional  supersonic  flows  with 
state-specific  vibrational  rate  kinetics.  These 
simulations  represent  some  of  the  energetics  ex¬ 
pected  to  occur  in  arcjets.  This  work  is  a  nec¬ 
essary  first  step  before  complete  modelling  of 
the  discharge  including  the  effects  of  electron- 
impact  processes  (see  also  (2]).  The  numerical 
method  used  here  is  the  Linearized  Block  Im¬ 
plicit  (LBI)  method  of  Briley  &  McDonald  {3]. 
The  LBI  method  has  excellent  stability  prop¬ 
erties  and  extends  easily  to  multi-dimensions. 
Formulation  of  the  boundary  and  initial  condi¬ 
tions  for  this  problem  are  discussed  in  detail, 
and  these  essentially  carry  over  for  similar  2-D 


and  3-D  cases. 


Nomenclat  ure 


A 

cross-sectional  area  of  the  nozzle  (m2) 

e 

internal  energy  per  unit  volume  (i/m3) 

Ei,v 

energy  of  vibrational  level  v 
of  species  i  (J) 

k 

Boltzmann’s  constant  (J/K) 

Mt.v 

molecular  mass  of  the  diatomic 
species  (kg) 

Mm 

atomic  or  molecular  mass  of  the 
diluent  species  (kg) 

n«,v 

number  density  of  vibrationally  active 
species  i  at  level  v  (m-3) 

number  density  of  diluent  species 
m  (m~3) 

Nd 

number  of  vibrationally  active  species 

Nm 

number  of  diluent  species 

P 

pressure  (N /m2) 

Po 

reservoir  pressure  (N/m2) 

t 

time  (s) 

T 

temperature  (K) 

To 

reservoir  temperature  (K) 

u 

velocity  (m/s) 

Kiwj r 

highest  vibrational  level  considered 
in  the  diatomic  species 

X 

axial  or  streamwise  coordinate  (m) 

P 

mass  density  (kg/m3) 

1  Introduction 

Arcjet  thrusters  are  electrothermal  devices 
where  a  propellant  stream  is  ohmically  heated 
and  subsequently  expanded  to  supersonic 
speeds.  These  thrusters  are  capable  of  robust 
operation  using  a  variety  of  propellants  includ¬ 
ing  monatomic  and  polyatomic  gases.  The 
energy  added  into  the  gas  is  distributed  into 
thermal  heating  (i.e.  increase  in  translational 
mode  temperature)  and  into  various  internal 
modes  (  electronic  excitation,  vibration,  and 
rotation).  However,  it  is  mainly  the  thermal 
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heating  which  is  subsequently  converted  into 
directed  kinetic  energy  (i.e.  useful  thrust)  in 
the  nozzle.  Thus,  any  energy  distributed  into 
the  internal  modes  which  is  subsequently  not 
recovered  in  translational  motion  is  lost  and  re¬ 
ferred  to  as  a  "frozen  flow  loss”.  This  work  is  a 
first  step  toward  predicting  frozen  flow  losses. 

The  internal  mode  energy  exchange  between 
molecules  in  a  discharge  not  only  results  in 
re-distribution  of  input  electrical  energy,  but 
also  determines,  moderates,  and  governs  elec¬ 
tron  energy  distributions.  Electrons  which  nor¬ 
mally  acquire  large  kinetic  energies  from  ap¬ 
plied  electric  fields  in  monatomic  gases,  are 
restricted  from  doing  so  in  plasmas  contain¬ 
ing  appreciable  amounts  of  molecular  species. 
This  is  because  inelastic  collisions  between  elec¬ 
trons  and  molecules  results  in  transfer  of  elec¬ 
tron  kinetic  energy  to  vibrational,  electronic,  or 
rotational  excitation.  Since  energy  exchange 
between  molecules  also  occurs  in  these  dis¬ 
charges,  the  well  known  non- Maxwellian  dis¬ 
tribution  such  as  Treanor  distribution[4]  occurs 
among  the  molecules  such  as  JV2,  thereby  lead¬ 
ing  to  highly  non- Maxwellian  electron  energy 
distributions.  Such  distributions  can  occur  at 
low  pressures  in  molecules  such  as  N2  and  CO 
(100  torr)  [5]  or  supra-atmospheric  pressures 
(20  atm.)  [6].  Consequently,  it  is  likely  that 
the  concept  of  an  electron  temperature  in  arc- 
jet  plasmas  is  meaningless  and  that  one  must 
resort  to  electron  energies  or  a  mean  electron 
energy  instead. 

Arcjet  thrusters  operating  at  low  or  moder¬ 
ate  powers  in  the  range  of  1  KW  to  30  KW 
on  ammon.a  (N /?3)  or  hydrazine  (H^N?)  are 
expected  to  produce  significant  amounts  of  N2 , 
NH,  and  H 2.  The  pressure  varies  from  100  torr 
in  the  nozzle  to  an  atmosphere  or  several  at¬ 
mospheres  in  the  constrictor  or  upstream  of 
the  constrictor.  Typically,  the  E/N  in  these 
thrusters  vary  over  a  wide  range:  10-23V  — 
m2  <  E/N  <  10~16V -m2  (or  10~19V -cm2  < 
E/N  <  10-,2V  -  cm2).  Molecular  discharges 


in  N2  have  been  studied  extensively  in  the  past 
(for  instance  see  (5))  over  a  range  that  is  par¬ 
tially  applicable  to  arcjet  flows.  Such  studies 
offer  valuable  insight  into  the  nature  of  the  elec¬ 
tron  energy  distribution  function.  For  instance, 
Fig.  1  from  [5]  shows  the  calculated  Electron 
Energy  Distribution  Functions  (EEDFs)  in  N2 
discharges  for  a  variety  of  E/N  ratios.  It  is  to  be 
emphasized  that  this  figure  is  a  semi-log  plot  so 
that  a  Maxwell-Boltzmann  distribution  would 
be  shown  as  a  sloping,  straight  line  on  such  a 
plot.  As  can  be  seen  from  Fig.l,  the  EEDFs  are 
markedly  non-Maxwellian.  This' has  dramatic 
consequences  for  evaluation  of  transport  prop¬ 
erties  since  cross-sections  for  collisional  pro¬ 
cesses  are  averaged  over  energy  distributions 
(usually  taken  to  be  Maxwellian). 

It  is  also  instructive  to  examine  the  "frac¬ 
tional  power  transfer”  versus  E/N  for  N2  dis¬ 
charges  as  given  by  Nighan  [5],  shown  here  in 
Fig.  2.  With  respect  to  the  arcjet,  this  figure 
may  be  interpreted  as  follows.  At  low  pow¬ 
ers  characterized  by  the  lower  E/N  values  dis¬ 
played  in  Fig.  2,  it  can  be  seen  that  much  of 
the  energy  added  to  the  discharge  goes  into  the 
vibrational  modes  of  N2.  Electronic  excitation 
gains  prominence  as  E/N  is  increased,  with  ion¬ 
ization  only  begining  to  consume  input  power 
at  E/Ns  approaching  10"15V  —  cm2.  Note  the 
conspicuous  absence  of  dissociation,  which  is 
understandable  since  the  dissociation  energy  of 
JV2  is  9.76  eV  and  the  E/Ns  in  Fig.  2  are  lower 
than  those  corresponding  to  average  electron 
energies  less  than  3  eV.  It  is  therefore  clear 
from  this  example  that  scaling  up  in  power  for 
arcjets  operating  on  nitrogen,  ammonia,  or  hy¬ 
drazine  is  strongly  dependent  on  the  degree 
of  dis-equilibria  between  the  various  internal 
modes.  Frozen  flow  losses  and  losses  in  effi¬ 
ciency  are  therefore  related  to  vibrational,  elec¬ 
tronic,  and  rotational  mode  non-equilibrium. 
Several  recent  works  have  already  pointed  to 
striking  differences  between  rotational,  vibra¬ 
tional,  and  translational  energies  in  the  plumes 
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of  arcjet  thrusters[l,  7]. 

Flows  with  high  degree  of  vibrational,  elec¬ 
tronic,  and  chemical  non-equilibrium  occur  in 
many  other  engineering  applications.  These 
include  flow  in  gas  dynamic  lasers,  rocket 
thrusters  (chemical  and  electrical)  and  chem¬ 
ical  vapor  deposition  (CVD)  processes.  A  typi¬ 
cal  feature  of  these  non-equilibrium  flows  is  the 
large  variation  in  characteristic  times  for  inter¬ 
nal  processes  which  can  vary  by  as  much  as 
thirty  to  forty  orders  of  magnitude  or  more  in 
various  regions  of  a  flow.  In  applications  where 
the  effects  of  state-specific  kinetics  or  large 
number  of  elementary  reactions  have  to  be  in¬ 
corporated,  there  is  a  need  for  a  robust,  effi¬ 
cient,  reliable,  and  accurate  numerical  method 
to  simulate  such  flows. 

Recently,  numerical  solutions  to  supersonic 
nozzle  flows  with  high  degree  of  vibrational 
non-equilibrium  and  state-specific  kinetics  have 
been  presented  [8].  These  quasi  1-D  solutions 
were  obtained  by  integrating  the  steady  forms 
of  the  governing  equations  using  the  stiff  equa¬ 
tion  integrator  LSODE  [9].  While  it  is  pos¬ 
sible  to  obtain  solutions  to  1-D  problems  us¬ 
ing  this  strategy,  it  does  not  extend  to  multi¬ 
dimensions.  This  is  a  serious  drawback  if 
one  wants  to  study  such  flows  in  two-  and 
three-dimensions,  while  including  the  effects  of 
heat  conduction  and  viscosity.  Our  interest  in 
these  multi-dimensional  flows  arises  from  ihe 
need  to  quantify  the  frozen  flow  losses  and  re¬ 
distribution  of  input  electrical  energy  in  arcjet 
thrusters.  Kence,  it  is  desirable  to  have  an  ef¬ 
ficient  method  that  can  handle  such  stiff  prob¬ 
lems  and  extend  from  1-D  all  the  way  to  3-D. 
The  Linearized  Block  Implicit  (LB I)  scheme  of 
Briley  &  McDonald  [3]  is  such  a  method.  As  we 
show  in  a  recent  work  [2],  this  method  is  very 
well  suited  for  studying  viscous,  internal  flows 
in  geometries  with  high  exit  to  throat  area  ra¬ 
tios.  Our  ultimate  goal  is  to  be  able  to  study 
these  internal  flows  after  including  the  effects  of 
vibrational  and  chemical  non-equilibrium  with 


state-specific  kinetics. 

In  this  paper,  we  tackle  the  first  steps  in 
calculating  vibrational  mode  non-equilibrium, 
simulating  the  energy  loading  in  the  arcjet 
thruster.  It  must  be  pointed  out  at  the  out¬ 
set  that  the  results  we  present  here  do  not 
account  for  the  electric  discharge  in  the  arc¬ 
jet  unlike  [2].  Rather,  we  simulate  the  energy 
loading  of  the  arc  as  an  initial  total  tempera¬ 
ture  to  which  the  propellant  stream  is  heated 
before  expansion  through  the  nozzle.  Work 
incorporating  vibrational/electronic  mode  dis¬ 
equilibrium  together  with  the  discharge  phe¬ 
nomena  is  presently  underway,  and  will  be  re¬ 
ported  later.  We  focus  here  on  flows  of  pure  N?, 
Ni/Hi,  Ni/ At,  and  Ni/CO/Ar  mixtures  to  il¬ 
lustrate  the  effects  of  inter- molecular  and  intra¬ 
molecular  V-V  (Vibration-to- Vibration)  and 
V-T  (Vibration-to- Translation)  energy  trans¬ 
fer.  Our  present  work  thus  differs  in  ap¬ 
proach  and  scope  from  previous  or  ongoing 
work  [10,  11,  12,  13,  14,  15,  16,  17]. 

This  paper  is  organized  as  follows.  Formu¬ 
lation  of  the  problem  which  differs  from  [8]  is 
discussed  in  the  next  section,  followed  by  a  dis¬ 
cussion  of  the  numerical  method  in  section  3. 
Results  are  discussed  in  section  4,  followed  fi¬ 
nally  by  a  summary  and  conclusions. 

2  Formulation 

In  the  arcjet,  much  of  the  ohmic  dissipation 
takes  place  in  the  constrictor  region  and  up¬ 
stream  of  the  nozzle.  This  addition  of  energy 
to  the  upstream  subsonic  propellant  can  be 
viewed  as  an  addition  of  thermal  energy  to  the 
flowing  propellant  stream  driving  it  to  choke 
and  accelerate  to  supersonic  speeds.  Although 
the  following  is  not  a  completely  true  descrip¬ 
tion  of  the  arcjet  operation,  the  present  results 
serve  to  illustrate  the  role  of  molecular  energy 
transfer  in  arcjet  flows.  Since  the  problem  for¬ 
mulation  and  solution  methodology  that  is  used 
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in  the  present  work  is  completely  different  from 
that  used  in  [8],  we  present  them  in  detail  here. 
The  governing  equations  are  the  quasi  1-D  gas 
dynamic  equations  together  with  spec '<  us  con¬ 
servation  equations  [8].  These  are  written  as 
follows  : 


9p  ,  1  djpuA) 
dt  A  dx  ’ 

d{pv)  1  d(pu7 A)  _  dp 
dt  +  A  dx  dx 


(1) 

(2) 


optically  thin  gases  considered  here.  These  as¬ 
sumptions  do  not  indicate  a  loss  of  generality, 
since  all  these  terms  can  be  included  in  the  for¬ 
mulation  if  applicable. 

The  pressure  is  related  to  the  temperature 
and  the  number  densities  through  the  following 
equation  of  state  : 

(Nd  Vm„  Nm  \ 

E  +  E n™ ) kT  •  (6) 

1=1  v=0  m=l  / 


d_ 

dt 


e  +  P  + 


=  0 
(3) 


dnitV  1  g(n,,„Ui4) 
dt  +A  dx 


VTitV+VVi<v+SRDi,v  , 


and  the  mass  density  is  related  to  the  number 
densities  as  follows  : 


^  =  E  E  +  E  .  (?) 

1=1  v=0  m= 1 


1  <  i  <  Nd,  0<v<  Vmax  (4) 


dnm  ,  1  d(nmuA)  n 
dt  +A  dx  ~ 


1  <  m  <  Nm 


(5) 

where  the  subscripts  i,t;  refer  to  diatomic 
species  *  in  vibrational  state  v  and  the  sub¬ 
script  m  refers  to  diluent  species  m  (usually 
monatomic).  In  addition,  only  the  diatomic 
species  are  assumed  to  be  vibrationally  active 
and  Vmax  +  1  is  the  total  number  of  vibra¬ 
tional  levels  considered  in  each  of  the  diatomic 
species.  is  calculated  by  treating  the  di¬ 
atomic  molecules  as  anharmonic  oscillators  [8]. 
The  heat  of  formation  of  various  species  can 
be  included  in  the  energy  equation.  However, 
since  chemical  reactions  are  not  considered  in 
the  examples  presented  here,  they  have  been 
left  out  of  Eq.  (3).  Also,  energy  of  the  elec¬ 
tronic  levels  is  neglected  as  we  do  not  include 
electronic  transitions  in  the  examples.  The  ra¬ 
diative  loss  which  would  appear  as  a  sink  term 
on  the  right  hand  side  of  Eq.  (3),  is  neglected 
here  since  it  is  expected  to  be  small  for  the 


where  M  stands  for  the  atomic  or  molecular 
mass  of  the  respective  species.  The  internal  en¬ 
ergy  per  unit  volume  is  given  by  the  following  , 
assuming  equilibrium  between  the  translational 
and  rotational  modes  : 


1  Mn  N d  Vmax  c 

m=l  1=1  v=0 

(8) 

The  terms  on  the  right  hand  side  of 
Eq.  (4)  represent  the  state-specific  Vibration- 
Translation  (VT),  Vibration- Vibration  (VV) 
and  Spontaneous  Radiative  Decay  (SRD) 
terms  as  described  by  Dunnwald  et  al.  [18]  and 
the  reader  is  referred  to  their  work  for  the  de¬ 
tailed  expressions  for  these  terms. 

2.1  Boundary  Conditions 

We  assume  that  a  large  reservoir  is  present  up¬ 
stream  of  the  nozzle  and  that  the  conditions 
in  the  reservoir  such  as  temperature,  pressure, 
and  volume  fractions  of  the  various  species 
present  are  known.  Further,  it  is  assumed  that 
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the  diatomic  species  are  in  Boltzmann  equilib¬ 
rium  both  in  the  reservoir  and  at  the  nozzle 
inlet  plane  and  that  the  mixture  is  composi- 
tionally  frozen  between  these  two  locations.  In 
addition,  heat  losses  and  frictional  losses  in  go¬ 
ing  from  the  reservoir  to  the  nozzle  inlet  are 
assumed  to  be  negligible,  without  any  loss  of 
generality. 

With  the  above  assumptions,  we  can  relate 
the  conditions  at  the  nozzle  inlet  to  those  at 
the  reservoir  by  conservation  of  energy.  This 
gives  us  an  equation  that  involves  the  veloc¬ 
ity,  temperature,  and  total  number  density  (= 
z!t\  Ew=o*  n*.v+£m= i  »m)  at  the  inlet.  Since 
the  volume  fractions  of  the  monatomic  species 
at  the  inlet  are  known,  we  can  relate  the  indi¬ 
vidual  nm’s  to  the  total  number  density.  For 
the  diatomic  species,  the  known  volume  frac¬ 
tions  at  the  inlet  together  with  the  fact  that 
they  are  in  Boltzmann  equilibrium  [19]  allows 
us  to  relate  the  »i>t,’s  to  the  total  number  den¬ 
sity.  Finally,  we  can  relate  the  inlet  static  pres¬ 
sure,  and  hence  the  total  number  density  (see 
Eq.  (6)),  and  temperature  to  their  respective 
values  at  the  reservoir  through  the  isentropic 
equation  of  state  [20].  In  the  present  case, 
the  velocity  at  the  inlet  of  the  nozzle  is  small 
because  of  the  large  inlet-to-throat  area  ratio 
and  thus  the  use  of  the  isentropic  equation  of 
state  is  acceptable.  These  conditions  allow  us 
to  eliminate  all  the  dependent  variables  except 
the  velocity  at  the  inlet,  which  is  determined 
through  an  implicit  extrapolation  [3].  At  the 
exit  to  the  nozzle,  implicit  extrapolation  is  used 
for  all  the  variables  since  we  are  interested  only 
in  supersonic  flows. 

Although  the  governing  equations  and 
boundary  conditions  have  been  presented  in 
their  dimensional  form  in  this  section,  it 
must  be  pointed  out  that  they  are  non- 
dimensionalized  before  being  solved  numeri¬ 
cally.  This  we  do  by  using  upstream  reservoir 
conditions  as  reference  quantities.  This  choice 
is  especially  convenient,  since  the  reservoir  con¬ 


ditions  need  to  be  specified  as  boundary  condi¬ 
tions  anyway  as  we  have  already  seen. 

3  Numerical  Method 

Equations  (l)-(5)  together  with  the  boundary 
conditions  are  discretized  in  time  and  linearized 
implicitly  using  the  procedure  outlined  by  Bri¬ 
ley  &  McDonald  [3].  Artificial  dissipation  in 
the  form  of  a  second  derivative  in  the  axial 
or  streamwise  coordinate  is  added  to  the  right 
hand  side  of  equations  (l)-(5)  to  ensure  numer¬ 
ical  stability  [3].  These  equations  are  then  dis¬ 
cretized  using  second  order  accurate  central  dif¬ 
ferences  on  a  non-uniform  grid  with  fine  spac- 
ings  near  the  throat  region.  After  eliminating 
the  boundary  points  from  the  resulting  system 
of  linear  algebraic  equations  using  the  bound¬ 
ary  conditions,  a  block  tridiagonal  system  of 
equations  results  which  can  be  solved  by  using 
standard  block  LU  decomposition  algorithms 
[21].  We  must  note  at  this  point  that  this  sys¬ 
tem  is  linearly  dependent  due  to  the  relation¬ 
ship  between  the  mass  density  p  and  the  num¬ 
ber  densities  through  Eq.  (7),  which  holds  at 
each  grid  point.  In  order  to  make  this  system 
linearly  independent,  we  arbitrarily  replace  one 
of  the  species  equations,  such  as  the  one  cor¬ 
responding  to  i=l  and  v~0  at  each  grid  point 
with  Eq.  (7)  applied  at  that  grid  point.  Note 
that  this  is  done  implicitly  so  that  overall  sta¬ 
bility  is  not  lost.  The  system  of  equations  is 
now  rendered  linearly  independent  and  ready 
to  be  solved. 

Since  we  are  solving  the  unsteady  form  of 
the  governing  equations  here,  the  frozen  flow 
solution  is  prescribed  as  the  initial  guess  to 
start  the  solution  process.  This  solution  is  ob¬ 
tained  by  numerically  solving  the  above  sys¬ 
tem  of  equations  with  the  right  hand  side  of 
Eqs.  (4)  set  identically  to  zero.  In  principle, 
a  truly  frozen  flow  cannot  be  used  as  an  initial 
guess  for  non-equilibrium  calculations.  How- 
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ever,  in  practice,  these  solutions  exhibit  devi¬ 
ations  from  the  truly  frozen  flow  solution  due 
to  discretization  and  round-off  errors  normally 
present  in  such  numerical  solutions.  Therefore, 
the  frozen  flow  solution  can  be  used  as  initial 
guess  for  the  non-equilibrium  calculations. 

4  Results  &  Discussion 

We  consider  the  geometry  of  the  nominal  30 
KW  arcjet.  However,  as  mentioned  earlier,  we 
do  not  model  the  arc.  Rather,  the  effect  of  the 
arc's  heating  of  the  flow  is  simulated  as  an  ef¬ 
fective  upstream  total  temperature.  The  cross- 
section  of  the  nozzle  as  a  function  of  the  axial 
coordinate  is  shown  in  Fig.  3.  For  all  the  re¬ 
sults  presented  here,  a  total  of  120  grid  points 
in  the  axial  direction  is  used.  The  unsteady, 
non-equilibrium  calculations  are  stopped  when 
the  percent  change  in  the  variables  is  less  than 
or  equal  to  0.00001  %.  The  size  of  the  dimen¬ 
sionless  time  step  is  about  5  x  10~7  seconds  for 
all  the  cases,  even  though  larger  values  can  be 
used  for  some  of  the  cases.  Since  the  number 
of  grid  points  is  the  same  for  all  the  cases,  the 
total  computing  time  for  each  case  depends  al¬ 
most  entirely  upon  the  the  size  of  the  blocks  in 
the  block  tridiagonal  system,  which  is  equal  to 
3  +  Nm  +  NdVmax-  In  the  present  study,  block 
sizes  as  large  as  84  (corresponding  to  Nm= 1, 
Nd= 2  and  Vrmai=40)  are  considered.  The  to¬ 
tal  computing  time  (including  the  time  taken 
for  generating  the  initial  guess)  for  this  case 
of  largest  block  size  is  about  4  minutes  in  uni¬ 
processor  mode  on  the  Cray  Y-MP,  while  it  is 
about  8  minutes  for  the  cases  with  a  single  di¬ 
atomic  species. 

We  first  consider  the  flow  of  pure  iV2  with 
the  reservoir  temperature  at  5000  K  and  the 
reservoir  pressure  at  2  atm.  Only  the  ground 
electronic  state  is  modeled.  The  mass  density, 
velocity  and  temperature  along  the  length  of 
the  nozzle  are  plotted  in  Fig.  4.  As  can  be 


seen  in  this  figure,  the  gradients  of  these  quan¬ 
tities  are  quite  large  near  the  constrictor,  and 
the  open  symbols  (shown  only  for  the  veloc¬ 
ity  for  the  sake  of  clarity)  clearly  show  that 
the  non-uniform  grid  is  able  to  resolve  these 
gradients  very  well.  Fig.  5  shows  normalized 
population  distributions  for  N?  at  three  sec¬ 
tions,  namely,  inlet,  throat  (exit  to  the  con¬ 
strictor)  and  the  exit.  The  distribution  at  the 
throat  is  very  close  to  a  Boltzmann  distribu¬ 
tion,  while  the  distribution  at  the  exit  shows 
considerable  departure  from  equilibrium.  As  is 
evident,  the  initial  energy  loading  of  the  gas  up¬ 
stream  (i.e.  via  a  high  stagnation  temperature) 
is  re-distributed  into  the  upper  vibrational  lev¬ 
els  as  the  gas  suddenly  expands  and  cools.  This 
effect  is  well  known,  and  referred  to  as  an- 
harmonic  V-V  (Vibration- Vibration)  pumping 
[4].  This  effect  occurs  when  conditions  are  such 
that  V-V  rates  are  much  faster  than  V-T  rates, 
typical  of  conditions  in  supersonic  expansions. 
These  distributions  agree  well  with  those  pre¬ 
sented  in  [8]. 

The  effects  of  adding  diluents  in  the  form 
of  monatomic  species  was  also  studied.  These 
diluents  are  added,  depending  on  the  applica¬ 
tion,  in  order  to  either  inhibit  or  promote  V- 
T  relaxation  downstream  of  the  constrictor  re¬ 
gion.  For  the  20%  JV2/ 80%  Ar  mixture  shown 
in  Fig.  6,  the  JV2  is  seen  to  be  less  V-V  pumped 
than  for  the  pure  N 2  case  (Fig.  5).  Evidently, 
the  vibrational  modes  of  Ni  are  de-excited  via 
V-T  transfer  to  argon.  Addition  of  10%  H i 
also  results  in  de-excitation  of  the  vibrational 
modes  of  JV2  comparable  to  the  20%1V2  /80% 
Ar  case,  as  can  be  seen  from  Fig.  7.  The  popu¬ 
lation  distributions,  however,  are  not  as  sensi¬ 
tive  to  small  changes  in  the  percentages  of  Hi 
in  the  mixture  as  was  reported  in  [8].  This  is 
due  to  an  error  in  the  values  of  the  constants 
Ahk,  Bj'k  and  Cj,k  (see  Table  1  in  de  Roany 
et  al.  [8])  used  for  calculating  the  CO- Hi  V-T 
rates.  The  correct  values  for  these  constants 
are  shown  in  Table  1. 
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As  an  illustration  of  inter-molecular  energy  brational  populations  of  the  ground  electronic 
transfer,  we  consider  a  mixture  of  20%CO,  states  of  the  relevant  diatomic  molecules.  Since 
20%  Wj,  and  60%  Ar.  In  this  mixture,  both  electron  energies  are  primarily  dependent  on 
CO  and  N2  are  vibrationally  active,  but  only  inelastic  collisions  with  molecules  such  as  N 2, 
their  ground  electronic  manifolds.are  modeled,  their  distribution  in  non-Boltzmann  or  non- 
We  neglect  vibrational  energy  transfer  to  low-  Maxwellian  since  these  diatomic  species  ex- 
lying  electronic  levels,  although  the  latter  could  hibit  significant  departures  from  equilibrium, 
also  be  included  within  the  framework  of  the  Hence  transport  properties  which  depend  on 
present  method.  The  vibrational  populations  cross  sections  averaged  over  these  distributions 
of  these  species  at  the  exit  plane  are  plotted  are  markedly  affected.  This  work  underscores 
in  Fig.  8.  The  distributions  show  consid-  the  importance  of  considering  inelastic  electron 
erable  non-equilibrium  and  the  number  den-  impact  processes  for  evaluating  transport  prop- 
sities  of  N2  for  vibrational  levels  above  5  erties  in  arcjet  flows  correctly, 
are  about  2  orders  of  magnitude  lower  than  These  solutions  are  obtained  with  the  LBI 
those  of  CO.  These  distributions  agree  quali-  scheme  of  Briley  &  McDonald  [3]  using  state- 
tatively  with  those  presented  in  [8].  The  aver-  specific  kinetics.  This  algorithm  has  excellent 
age  vibrational  energy  per  species,  defined  as  stability  properties  and  allows  the  use  of  quite 
53«=o*  ni,vEi,v  /  n«.«  's  Pitted  in  Fig.  large  time  steps  even  in  such  stiff  systems. 

9  for  CO  and  N2.  Dramatic  transfer  of  vibra-  This  feature  is  particularly  desirable  when  per¬ 
sonal  energy  from  N2  which  has  more  widely  forming  such  calculations  in  multi-dimensions, 
spaced  energy  levels,  to  CO  which  has  more  The  boundary  and  initial  conditions  that  are 
closely  spaced  levels  can  be  seen  in  this  figure,  discussed  here  carry  over  to  multi-dimensional 
This  is  in  agreement  with  theoretical  predic-  problems  as  well.  The  present  method  is  there- 
tions  [4]  as  well  as  the  numerical  results  pre-  fore  well  suited  to  studying  high  speed,  multi- 
sented  in  [8].  When  ammonia  or  hydrazine  dimensional,  reacting,  plasma  flows  such  as 
is  used,  species  such  as  N2,  H2,  and  NH  are  those  encountered  in  arcjet  thrusters,  lncor- 
produced  in  the  discharge.  Thus  vibrational  poration  of  the  detailed  knietic  processes  dis- 
energy  transfer  can  take  place  between  these  cussed  herein  into  a  2-D  axi- symmetric  model 
species,  ultimately  affecting  the  electron  en-  reported  in  [2]  is  presently  underway, 
ergy  distribution  and  making  it  markedly  non- 
Maxwellian.  The  presence  of  such  vibrationally 
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the  inlet  stagnation  temperature.  The  results 
are  therefore  somewhat  representative  of  the  vi- 
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APPENDIX 


The  detailed  expressions  for  the  VT,  VV  and 
the  SRD  terms  are  given  here.  Even  though 
the  expressions  presented  here  look  identical  to 
those  given  in  [8],  there  are  some  minor  dif¬ 
ferences  between  the  two  due  to  typographical 
errors  in  [8]. 

Ei,v,  the  vibrational  energy  at  level  v  of 
species  i  is  given  by  : 


■Ei.v  — 


1.6022  x  10~19 
8065.479 


Also,  Pjj1  1  is  the  rate  constant  (in  m3/s)  of 
the  VT  transition  : 


where  Xi(t>)  refers  to  species  t  in  vibrational 
state  v.  The  rate  constant  is  given  by  [8], 


y,Jl-Xeiv ’ ’ 


where 


FW=U 3-e~2^3)  e"2^3 


+  D  j 

where  the  electronic  energy  has  been  ne¬ 
glected  since  we  consider  only  ground  electronic 
states  here.  In  general,  is  also  depen¬ 
dent  on  the  rotational  quantum  number  [22], 
but  in  this  work  we  treat  the  rotational  modes 
as  being  in  equilibrium  with  the  translational 
modes.  (in  cm-1)  and  \ei  (dimensionless) 
are  the  spectroscopic  constants  for  molecules  of 
species  i  (see  Table  4). 

Vibration-Translation  (VT)  term 

The  VT  term  in  the  species  equation  (4)  can 
be  written  in  form  as 

Nd+Nm 

VT...  = 

j- 1 

-exp  {-&E?/kT)ni,v) 

-p?r'  («i,,-«p(-A£r-,/trK„-1)] 


with 


_  0-3/2 
AiJ  ~  i 


,/Sdagq 

V  t  *e. 


Here  A E"_1  =  E,tV  -  E.^-i  is  the  differ¬ 
ence  in  energy  between  products  and  reactants 
in  the  VT  transition;  0,  =  hcu>ei/k  (with  c, 
the  speed  of  light  in  cm/s)  is  the  vibrational 
characteristic  temperature  (in  K)  of  species  i; 
Qij  =  16xVi .jC^ei^/k  is  in  K  where  is 
the  reduced  mass  in  kg  and  /  (=  0.2  x  10_1°  m) 
is  the  range  parameter  [23]. 

In  the  above  expression  for  the  rate  constant, 
Vij  is  a  coefficient  that  allows  for  fitting  to 
the  available  experimental  relaxation  data  for 
the  1-0  transition  rate  i.e  .P’*’0.  It  is  given  (in 
m3/s)  as  : 


10~13(1  ~  XeiW 


Here  (rp)t  J  is  the  vibrational  relaxation  time 
in  ns.atm  and  is  expressed  as  : 


where 

Vmar 

=  J2  »>.« 

ui=0 

in  the  case  of  a  vibrationally  active  species 
and  nj  =  nm  in  the  case  of  a  diluent  species. 


+Cij(T)~V3 

where  Aij ,  Bij  and  Qj  are  empirical  con¬ 
stants  chosen  to  match  the  experimental  data. 
They  are  given  in  Table  1 . 
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Vibration- Vibration  (W)  term 


The  VV  term  in  the  species  equation  (4)  can 
be  written  in  as 


vk> = e  E  ci:—  (».>«»,> 

;=1  wsO 


-  exp  (-  (A£V  -  A£”)  /fcr)  »<.*»*„+!  ] 


_Ov*v-i 

ij;u),ti)+l 


where  (in  A'-1),  £tJ  in  K  and  6fJ  in 
K  are  empirical  parameters  fitted  to  available 
experimental  data.  Their  values  are  given  in 
Table  2.  Also,  F( A)  is  the  same  as  before  in 
the  case  of  VT  transfer  while  the  expression 
for  A  now  becomes 


«u,v-l 

1,id 


2-3/2 


[e~\AEri-*Erl\ 
V  t  ke. 


Finally 


-exp  (-  (A ETl  ~  ±E?)  /&)  ni.„_1nitU,+1] 

where  is  the  rate  constant  (in 

m3/s)  of  the  V-V  transition  : 


Xi(v)  +  Xj(w  -  1)  Xi{v  -  1)  +  Xj( w). 


The  rate  constant  is  written  as  a  sum  of  a 
short  range  and  a  long  range  contribution  [8] 
as  follows  : 

QV,V-1  _  n  / 

l.ui  "*J  l,ui 


.  T  \  .-( AB,—1  - )/2kT 

I’t'iJ;  w~l,w)e  3 

with 

=  IQ"6  x  400a?jy/*kT/2»ij 


where  ZiJ  is  the  collision  number  in  m3/s 
and  Tcafi  is  the  collision  cross-section  in  cm2 

*W  _ 

with  o  =  3.75  x  10  8  cm.  We  also  have 


and 


c«.v— 1  _  r  t  *•’ 

~  l~X~V 


W 


1  -  Xrjw 


rti.V- 1 

l,w 


a,  -I-  1  N2 
+  3-2u/ 


t?(o,  +  2  -  2t>)(a,  +  4  -  2t>) 
a, (a,  +  3  -  v) 
where  a,  =  l/*e,. 


Spontaneous  Radiative  Decay  term 


The  SRD  term  in  the  species  equation  (4)  for 
radiative  transition  between  levels  of  the  same 
electronic  species  can  be  written  as  [8] 


S  RD^ 


=  I>: 

tl/=l 


t/+tu,v 


**t,V+ W  -4,’ 


v,v— w 


n, 


ttv 


where,  A ,  the  Einstein  A  coefficients  are 
given  by 


Avj'v~w  1 

A1’0  ~  (a«  -2)(at-3)U; 

tti  bv(bv  +  w)(bv  +  2w) 
(w  -  w)\  n ^(m-v  +  q) 


with  bv  =  ai  -  2v  -  1  and  w1  and  A)'0  (in 
1/s)  are  given  in  Table  3. 
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TABLE  1. 


VT 


Species 

Aij 

B-j  (A'1") 

C-j  (K3/3) 

ejye*2  (tf-1) 

Ref. 

CO-CO 

-15.23 

280.5 

-549.6 

45.6  x  10~2 

[24] 

CO- Ar 

10.38 

0.0 

0.0 

53.72  x  10~2 

[24] 

CO -h2 

-26.22 

-10.61 

-12.80 

6.089  x  10~2 

[25] 

Nr  Hi 

-26.22 

-10.61 

-12.80 

6.089  x  10~2 

[25] 

con2 

-7.934 

147.7 

0.0 

45.6  x  10~2 

[8] 

N2-N2 

-12.539 

258.9 

-390.9 

45.6  x  10“2 

[8] 

N2-At 

-15.62 

168.95 

0.0 

53.72  x  10“2 

[8] 

TABLE  2. 


VV 


Species 

w 

b-j(K) 

Ref. 

CO-CO 

1.64  x  10"6 

1.6142 

40.36 

[23] 

CON2 

7.006  x  10-8 

1.897  x  10-2 

191.42 

[8] 

N2-N2 

9.37  x  10-8 

0.0 

- 

[8] 

TABLE  3. 


SRD 

Species  w'  A*1,0  (1/s)  Ref. 

CO-CO  4  30.3  [23] 

N2-N2  -  "  0.0  [8] 

TABLE  4. 

Spectroscopic  Constants 

Species  w«X«  ©’  (A")  Ref. 

(cm-1)  (cm-1) 

CO-CO  2169.8  13.288  3121.0  [22; 

N2-N2  2358.6  14.324  3396.8  [22' 
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Fig.  9:  The  average  vibrational  energy  per  particle  (normalized  by  6.903  x  10  P 
is  plotted  here  versus  stream  wise  coordinate  (normalized  by  0.01  m)  for 
flow  of  N2/CO/Ar.  Note  how  the  nitrogen  up-pumps  the  CO  via 

vibrational  energy  transfer  and  is  de-exdted  in  the  process. 
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ABSTRACT 


Swirl  has  long  been  thought  to  be  necessary 
to  sustain  arcs  in  high-speed  flows.  In  this 
paper,  we  explore  the  effect  of  injected  swirl 
on  both  the  cold  flow  (i.e.  non-ionized  and 
non-reacting)  and  the  plasma  flow  in  arcjet 
thrusters.  The  fully  implicit  Linearized  Block 
Implicit  (LBI)  scheme  of  Briley  k.  McDonald[l] 
is  used  to  develop  highly  accurate  solutions  to 
internal,  viscous,  supersonic  flows  of  the  type 
encountered  in  arcjet  thrusters.  In  the  case  of 
cold  flow,  it  is  found  that  the  injected  swirl  per¬ 
sists  downstream  for  low  mass  flow  rates  (on 
the  order  of  mg/s)  and  is  actually  enhanced  at 
the  constrictor  entrance,  near  the  cathode  tip. 
However,  for  higher  mass  flow  rates  (on  the  or¬ 
der  of  g/s),  the  injected  swirl  is  overcome  by 
the  torque  due  to  viscous  forces  (both  in  the 
flow  as  well  as  at  the  walls)  and  fails  to  persist 
far  enough  downstream.  Consequently,  swirl 
stabilization  of  the  flow  in  arcjets  can  only  suc¬ 
ceed  under  certain  operating  conditions.  The 
actual  plasma  flow  is  also  investigated  in  the 
presence  of  swirl,  but  our  present  simulations 
show  continuous  decay  of  injected  swirl.  As 
we  show  in  this  paper,  this  is  an  artefact  and 
limitation  of  our  present  simulation.  Real  ar¬ 
cjet  flows  should  experience  the  effects  we  ob¬ 
serve  in  our  cold  flow  simulations.  This  effect 
has  important  ramifications  for  the  design  and 
locations  of  the  swirl  injectors  and  overall  ar¬ 
cjet  performance.  On-going  improvements  to 
our  simulation  to  rectify  this  problem  are  dis¬ 
cussed. 

Nomenclature 

Be  magnetic  induction  (T) 

e  internal  energy  per  unit  volume  (J/m3) 

jT  radial  component  of  current 
density  (A/m2) 

jx  axial  component  of  current 
density  (A/m2) 

k  thermal  conductivity  (W/m2  K) 


ka  Boltzmann’s  constant  (J/K) 

m0  mass  of  an  argon  atom  (kg) 

na  number  density  of  atoms 

(particle  /  m3) 

ne  number  density  of  electrons 

(particle  /  m3) 
p  pressure  (N/m2) 

r  radial  coordinate  (m) 

t  time  (s) 

T  temperature  (K) 

u  radial  velocity  (m/s) 

t;  swirl  velocity  (m/s) 

w  axial  velocity  (m/s) 

2  axial  coordinate  (m) 

€,  ionization  potential  (J) 

p  coefficient  of  viscosity  (Ns/m2) 

p0  permeability  of  free  space  (Tm/A) 
p  mass  density  (kg/m3) 

rfr  magnetic  stream  function  (Tm) 

a  electrical  conductivity  (mho/m) 


1  Introduction 

Despite  the  near  flight-readiness  of  arcjet 
thrusters,  the  plasma  flow  in  these  devices 
is  only  begining  to  be  understood.  Un¬ 
like  its  distant  cousin,  the  Magnetoplasma- 
dynamic  (MPD)  thruster,  the  arcjet  is  basi¬ 
cally  an  electrothermal  device.  The  propel¬ 
lant  (either  monatomic  or  molecular)  stream 
is  heated  by  an  electric  arc  and  expanded 
though  a  converging-diverging  nozzle  to  super¬ 
sonic  speeds.  Although  the  device  is  simple  in 
structure,  the  phenomena  that  govern  flow  in 
the  device  are  quite  complicated  and  are  yet 
to  deciphered.  Among  these  is  the  question  of 
whether  or  not  injected  swirl  is  necessary  for  ef¬ 
fective  operation,  and  whether  this  swirl  plays 
any  role  in  arcjet  fluid  dynamics  in  the  first 
place.  Swirl  in  the  flow-field  forms  the  central 
focus  of  this  paper. 

The  compressible  Navier-Stokes  equations 
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govern  the  physics  of  flow  fields  in  arcjets. 
These  flows  are  often  times  further  compli¬ 
cated  by  the  effects  of  departure  from  ther¬ 
modynamic  equilibrium  of  the  internal  modes 
of  molecular  motion.  This  disequilibrium  be¬ 
tween  internal  modes  (vibration,  rotation  and 
electronic  excitation)  and  external  modes  i.e. 
translational  molecular  motions  gives  rise  to 
complicated  chemistry.  Many  numerical  meth¬ 
ods  have  been  devised  to  solve  the  compress¬ 
ible,  Navier-Stokes  equations,  the  most  popular 
of  these  being  the  MacCormack  scheme  [2],  the 
Beam- Warming  scheme  [3]  and  the  Linearized 
Block  Implicit  (LBI)  scheme  of  Briley  &  Mc¬ 
Donald  [1].  Our  ultimate  aim  is  to  study  arcjet 
flows  including  the  effects  of  finite-rate  chem¬ 
istry  and  internal  mode  disequilibrium.  We 
have  chosen  the  LBI  method  for  this  purpose. 
This  method  has  the  advantages  of  being  able 
to  handle  large  numbers  of  species  as  well  as 
being  extendable  from  quasi-ID  through  2-D, 
to  fully  3-D  with  relative  ease. 

Recently,  several  groups  have  begun  numeri¬ 
cal  simulations  of  arcjet  flows[4, 5, 6,  7, 8, 9, 10]. 
Perhaps  the  earliest  of  these  simulations  were 
those  of  Keefer  et.  al.[4]  and  Butler  et.  al.[6,  7]. 
Both  groups  have  identified  o,  the  electri¬ 
cal  conductivity,  as  the  single  most  important 
property  that  governs  arcjet  flow.  This  is  be¬ 
cause  a  varies  by  orders  of  magnitude  between 
regions  in  the  arc  and  regions  near  the  colder 
electrode  boundaries,  a  determines  the  current 
distribution  in  the  gas,  and  hence  the  ohmic 
heating.  The  Ohmic  heating  in  turn  determines 
temperatures  and  ionization  fraction.  The  lat¬ 
ter  has  a  profound  influence  on  other  proper¬ 
ties,  most  notably  the  viscosityfll]  so  that  the 
device  performance  is  inevitably  altered.  Fur¬ 
thermore,  the  distribution  of  a  profoundly  af¬ 
fects  the  heat  addition  not  only  in  the  subsonic 
portion  of  the  flow  but  also  in  the  supersonic 
portion.  Although  the  latter  is  unavoidable, 
it  is  important  to  minimize  heating  the  super¬ 
sonic  flow  as  specific  impulse  will  suffer.  In 


addition,  no  existing  work  has  systematically 
examined  the  effects  of  injected  swirl  on  the 
flow-field  or  conductivity  distribution  in  arc- 
jets. 

In  this  paper,  we  systematically  examine  the 
effects  of  injected  swirl  using  an  axi-symrretric 
2-D  model  of  a  nominal  30  KW  arcjet.  Both 
cold-flow  and  flow  including  the  arc  are  exam¬ 
ined.  We  find  that  the  key  to  obtaining  so¬ 
lutions  to  these  flows  is  the  proper  modelling 
of  the  boundary  conditions  and  specification  of 
the  proper  initial  conditions  so  that  the  flow 
can  accelerate  smoothly  from  subsonic  to  su¬ 
personic  speeds. 

This  paper  is  organized  as  follows.  The 
governing  equations  and  boundary  conditions 
are  presented  first,  followed  by  the  numeri¬ 
cal  method.  As  brief  illustration  of  code  ver¬ 
ification,  sample  solutions  comparing  against 
experiments  solutions  are  presented  for  a 
converging-diverging  nozzle  similar  but  not 
identical  to  the  nominal  30  KW  arcjet  geom¬ 
etry.  Cold-flow  and  2-D  axisymmetric  arcjet 
solutions  with  and  without  swirl  are  discussed 
next,  followed  by  a  summary  and  conclusions 
of  this  work. 

2  Formulation 

The  equations  that  govern  the  cold  (i.e.  non¬ 
reacting)  viscous,  supersonic  flow  are  the  com¬ 
pressible  Navier-Stokes  equations.  For  the  arc¬ 
jet,  these  equations  become  considerably  more 
complicated  by  the  presence  of  a  body  force, 
heat  generation,  species  diffusion,  chemical  re¬ 
action,  ionization  and  recombination,  and  vari¬ 
able  properties.  However  in  this  paper,  we  con¬ 
cern  ourselves  only  with  a  monatomic  propel¬ 
lant  (Argon)  and  defer  consideration  of  molec¬ 
ular  propellants  in  a  companion  paper[12j. 
We  further  neglect  species  diffusion  terms,  as¬ 
sume  quasi-neutrality,  and  take  the  electron 
and  heavy-particle  temperatures  to  be  equal. 
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While  significant  departures  from  the  latter  are 
known  to  occur  in  the  colder  regions  of  the 
flow[ll],  we  are  mainly  interested  here  in  the 
flow  away  from  the  electrode- adjacent  regions. 
The  governing  equations  in  cylindrical  polar  co¬ 
ordinates  form  are: 
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and  t  is  the  internal  energy  per  unit  volume 
given  by: 


*  =  2^2”e  +  +  ”e€* 
and  the  p  is  the  pressure  given  by: 

p  =  (2ne  +  na)kgT 

^  In  the  above  equations, 
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where,  i4  =  ylrer  +  A$eg  +  Azez  is  any  vector 
and  <t>  is  a  scalar. 

In  addition,  we  have  the  following  species 
rate  and  magnetic  diffusion  equations: 
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(8) 


V  x  (V  x  vT/r) - x  (V  x  0/r) 

o 

-  n0(T  (v  x  (u  x  1 £/r))  =  0.  (9) 

The  thermal  conductivity,  coefficient  of  vis¬ 
cosity,  and  electrical  conductivity  are  com¬ 
puted  from  mean  free  path  theory  and  are  given 
by[13].  Here  kj  and  kr  are  the  forward  and  re¬ 
verse  rates  [11]  of  the  reaction 

Ar  +  e~  ^  Ar+  +  e~  +  t~  . 

The  mass  density  p  in  the  above  equations  is 
related  to  the  number  densities  as  follows  : 

p  =  ma{ne  +  na)  (10) 

where  ma  is  the  mass  of  a  neutral  atom. 

2.1  Boundary  Conditions 

Boundary  conditions  for  the  eight  dependent 
variables  ( p ,  u,  v,  w,  p ,  ip,  ne,  T)  need  to  be  spec¬ 
ified  along  the  inlet  and  exit  boundaries,  on  the 
wall,  and  along  the  centerline  of  the  nozzle  ge¬ 
ometry. 

No-slip  conditions  are  imposed  on  all  solid 
boundaries.  In  the  case  of  the  ionizing  arcjet 
flows,  the  temperature  at  the  anode  is  held  con¬ 
stant  at  1500  K,  while  the  cathode  temperature 
is  held  at  1000  K  along  its  straight  portion  and 
ramped  linearly  to  2000  K  along  its  sloping  fac 
up  to  the  tip.  The  mass  density  at  the  wa 
is  obtained  by  applying  the  radial  cornponen 
of  the  momentum  equation  on  the  wall.  The 
choice  of  the  radial  component  instead  of  the 
normal  component  is  dictated  by  the  need  for 
simplicity.  The  electron  number  density  on  the 
wall  is  obtained  by  setting  its  radial  derivative 
to  zero  there. 

Symmetry  conditions  are  imposed  along  the 
centerline  for  all  the  variables  except  the  swirl 
and  radial  components  of  velocity  v  and  u, 
which  are  set  equal  to  zero,  identically. 


Prescribing  boundary  conditions  at  the  in¬ 
let  and  exit  planes  of  the  nozzle  depends  on 
whether  the  flow  is  subsonic  or  supersonic.  All 
the  solutions  presented  here  consider  subsonic 
inlets.  Hence,  following  Briley  &  McDonald 
[1],  the  stagnation  temperature  and  stagnation 
pressure  at  the  inlet  are  specified.  In  addition, 
at  the  inlet,  the  radial  component  of  velocity  is 
set  to  zero  and  the  swirl  component  of  velocity 
is  either  zero  or  set  to  some  prescribed  distri¬ 
bution  corresponding  to  injection.  It  is  impor¬ 
tant  to  note  here  that  the  axial  component  of 
velocity  which  determines  the  mass  flow  rate 
through  the  nozzle,  is  not  specified  at  the  inlet 
but  is  determined  from  an  implicit  extrapola¬ 
tion  of  the  interior  values,  as  is  ne.  In  addition, 
the  mass  density  at  inlet  is  determined  by  using 
the  isentropic  relation  between  the  stagnation 
and  static  pressures  and  temperatures  there. 

Along  the  exit  plane,  the  flow  is  subsonic- 
supersonic.  All  the  variables  (including  ne)  ex¬ 
cept  mass  density  on  the  exit  plane  are  evalu¬ 
ated  by  implicit  extrapolation  from  the  interior 
[1].  If  the  local  Mach  number  A/defined  as  w/a , 
where  a2  =  |(1  +  ;J+no  )^bT /ma  is  the  frozen 
speed  of  sound  is  greater  than  one  at  a  point 
on  the  exit  plane,  then  the  mass  density  is  also 
evaluated  by  extrapolation.  If,  on  the  other 
hand,  M  <  1  at  some  point  on  the  exit  plane, 
then  the  exit  static  pressure  is  specified  and  the 
mass  density  is  evaluated  from  the  equation  of 
ate.  The  correct  value  of  the  exit  static  pres- 
e  would  usually  not  be  known,  unless  the 
ame  is  modelled  as  well.  Care  must  be  ex¬ 
ercised  in  choosing  a  value  for  the  exit  static 
pressure,  for,  if  it  is  too  high,  the  boundary 
layer  on  the  wall  in  the  diverging  section  of  the 
nozzle  can  separate.  This  is  an  important  detail 
that  often  escapes  mention  in  the  literature. 

The  boundary  conditions  on  ip  could  use  im¬ 
provement.  In  this  work,  we  prescribe  and  hold 
ip  fixed  at  a  certain  value  based  on  the  total  cur¬ 
rent  [4]  on  the  anode  wall  from  the  inlet  until 
the  exit  to  the  constrictor.  Beyond  this  point, 
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V>  is  varied  linearly  along  the  anode  wall  to  zero 
at  the  exit.  ^  is  taken  to  be  identically  zero  at 
the  exit  plane  and  at  the  centerline  (i.e.  r=0). 
Along  the  inlet  plane,  0  is  varied  linearly  from 
its  value  on  the  anode  wall  to  0.99  times  that 
values,  and  maintained  at  that  value  along  the 
straight  portion  of  the  cathode  as  well  as  on  the 
sloping  face  of  the  cathode  until  5  grid  points 
from  the  tip.  For  these  5  grid  points  near  the 
cathode  tip,  V'  is  decreased  linearly  from  0.99 
times  the  maximum  value  to  zero  in  order  to 
force  most  of  the  current  contours  to  approach 
the  cathode  axially. 

2.2  Initial  Conditions 

Due  to  the  wave  nature  of  the  governing  equa¬ 
tions,  initial  conditions  must  be  prescribed 
with  care.  In  general,  arbitrary  initial  values 
for  the  dependent  variables  will  not  yield  sta¬ 
ble  solutions  for  any  numerical  method.  Fur¬ 
thermore,  since  it  is  impossible  to  simulate  the 
breakdown  process  that  initiates  the  discharge 
ne  cannot  be  expected  to  start  from  0  (i.e.  com¬ 
pletely  non-ionized  state).  In  this  work,  we 
prescribe  an  upstream  stagnation  temperature 
(also  equal  to  the  reference  temperature  used 
for  non-dimensionalization)  of  10,000  K,  an  up¬ 
stream  stagnation  pressure  (which  depends  on 
the  desired  mass  flow  rate),  u  =  0,  v  —  0  (ex¬ 
cept  at  the  inlet,  if  swirl  is  injected),  w  =  0.1ao 
where  aa  is  the  frozen  speed  of  sound  evalu¬ 
ated  at  the  inlet  tagnation  temperature,  and 
ne  =  0.1  nrej  where  nrej  is  computed  from  the 
stagnation  pressure  and  stagnation  tempera¬ 
ture. 

It  is  evident  at  the  outset  that  an  unre¬ 
alistically  high  value  must  be  prescribed  for 
the  initial  stagnation  temperature  and  elec¬ 
tron  number  density.  This  is  a  drawback  of 
the  present  formulation.  It  suffices  to  state 
at  this  time  that  for  supersonic  flows  of  the 
type  encountered  in  arcjets,  the  upstream  in¬ 
let  static  pressure,  temperature,  and  densities 


must  adjust  themselves  in  order  to  accommo¬ 
date  the  presence  of  the  arc  downstream.  In 
other  words,  some  level  of  ionization  (and  tem¬ 
perature  comensurate  with  this  initial  value  of 
ionization)  must  be  prescribed  to  start  the  cal¬ 
culation.  Obviously,  the  final  converged  steady 
state  solutions  must  be  independent  of  the  ini¬ 
tial  value  chosen  for  ne.  In  this  work,  we  use  a 
simple  way  to  calculate  the  static  temeprature 
at  the  inlet  from  the  stagnation  temperature. 
As  a  result  of  this  simplicity,  our  inlet  tem¬ 
peratures  must  be  necessarily  high  in  order  to 
conserve  overall  energy  in  the  converged  state. 
Physically,  this  alludes  to  the  fact  that  some 
concentrated  region  of  Ohmic  heating  must  ex¬ 
ist  upstream  of  the  constrictor  that  pre-ionizes 
the  gas  before  it  enters  the  main  arc  region. 


3  Numerical  Method 

Equations  (l)-(5),(8)  are  solved  a  non- 
staggered  non-uniform  grid.  This  grid  in  the 
physical  domain  is  transformed  into  a  uniform 
grid  in  a  computational  domain  using  an  ana¬ 
lytical  coordinate  transformation.  The  trans¬ 
formation  is  such  that  the  nozzle  contour  is 
transformed  into  a  rectangle  in  the  computa¬ 
tional  domain  (14,  section  5-6.1].  The  govern¬ 
ing  equations  together  with  boundary  condi¬ 
tions  are  transformed  to  the  computational  do¬ 
main  and  solved  there. 

The  transformed  forms  of  the  governing 
equations  (except  for  the  magnetic  diffusion 
equation)  are  solved  using  the  LBI  scheme  of 
Briley  &  McDonald  [1].  Since  all  the  details 
of  the  scheme  have  already  been  reported  else¬ 
where  [1],  they  willl  not  be  given  here.  At  each 
time  step  in  the  time-marching  sequence,  the 
magnetic  diffusion  equation  is  solved  using  odd- 
even  point  SOR  with  the  relaxation  factor  w 
set  to  1.0.  Thus,  the  ^  field  is  solved  for  in¬ 
dependently  at  each  time  step.  The  other  cou¬ 
pled,  linear  algebraic  equations  resulting  from 
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the  other  governing  equations  are  then  solved 
by  using  the  Douglas-Gunn  ADI  method.  The 
ADI  method  itself  involves  the  solution  of  block 
tridiagonal  systems  in  each  coordinate  direc¬ 
tion,  with  block  size  equal  to  the  number  of 
dependent  variables,  which,  in  our  case  is  eight. 
These  systems  can  be  solved  using  existing  effi¬ 
cient  LU  decomposition  methods  [15].  Follow¬ 
ing  Briley  &  McDonald  [1],  we  add  dissipation 
in  the  form  a  second  derivative  in  the  axial  co¬ 
ordinate  to  the  unsteady  governing  equations. 

4  Cold  Flow 

4.1  C-D  Nozzle 

In  this  section,  we  present  viscous  2-D  axi- 
symmetric  solutions  for  flow  in  a  C-D  noz¬ 
zle.  These  are  for  cold-flow  (i.e.  non-reacting 
and  non-ionizing)  of  air  (7  =  1.4,  n  =  18  x 
10"6  Ns/m2,Pr  =  0.701)  and  are  presented 
for  the  purposes  of  code  verification.  The  noz¬ 
zle  geometry  is  identical  to  the  geometry  of 
the  30  kW  Arcjet  thruster  (shown  in  Fig.  1) 
but  without  the  inlet  plenum  and  the  constric¬ 
tor.  The  transition  from  the  converging  to  the 
diverging  section  is  achieved  within  five  grid 
points  through  the  use  of  a  quadratically  vary¬ 
ing  profile  for  the  outer  radius.  The  exit  to 
throat  area  ratio  for  this  nozzle  geometry  is 
25.  Solutions  presented  here  are  obtained  on 
a  non-uniform  grid  with  256  points  in  the  axial 
direction  and  64  points  in  the  radial  direction. 
The  most  interesting  feature  of  this  flow  is  the 
presence  of  weak  compression  waves  in  the  di¬ 
verging  section.  Profiles  of  Mach  number  along 
the  radial  grid  lines  j  =  0  and  j  =  38  are  shown 
in  Fig.  2  together  with  the  Mach  number  pro¬ 
file  from  the  corresponding  frictionless,  quasi 
1-D  solution.  Here,  the  grid  line  corresponding 
to  j  =  0  represents  the  centerline.  The  profile 
along  the  j  =  38  grid  line  is  closest  to  the  quasi 
1-D  profile.  The  maximum  exit  Mach  number 
in  the  viscous  case  is  about  4.6,  which  is  less 


than  the  value  of  about  5.0  for  the  frictionless, 
quasi  1-D  case.  The  profile  along  the  centerline 
shows  a  sharp  drop  in  the  Mach  number  near 
a  dimensionless  axial  coordinate  of  8.  One  can 
see  from  the  Mach  number  contours  in  Fig.  3,  a 
change  in  the  shapes  of  these  contours  near  the 
same  axial  location  (z  =  8).  In  this  figure,  the 
sonic  line  (M=  1)  is  shown  as  the  thick  line.  The 
radial  Mach  number  profiles  have  their  maxima 
located  on  the  centerline  up  until  z  =  8  and 
these  maxima  move  to  a  location  near  the  cen¬ 
ter  of  the  channel  once  we  cross  this  particular 
axial  location.  These  contours  look  very  similar 
to  those  obtained  using  the  method  of  charac¬ 
teristics  [16]  for  inviscid  flow  in  a  similar  ge¬ 
ometry.  They  also  agree  very  well  qualitatively 
with  the  experimental  results  presented  in  [17]. 
It  was  first  shown  by  Darwell  &  Badham  [16] 
that  compression  waves  could  be  generated  in 
C-D  nozzle  geometries  from  the  location  where 
the  curved  throat  section  meets  the  diverging 
conical  section,  if  the  radius  of  curvature  of  the 
nozzle  wall  profile  is  discontinuous  at  this  loca¬ 
tion.  These  compression  waves  can  coalesce  at 
the  centerline  at  a  point  downstream  to  form 
a  weak  oblique  shock.  Experiments  conducted 
subsequently  by  Back  et  al.  [18,  19]  and  Cuf- 
fel  et  al.  [17]  showed  this  to  be  indeed  true. 
The  sharp  drop  in  the  Mach  number  along  the 
centerline  in  Fig.  2  is  due  the  oblique  shocks  in¬ 
tersecting  the  centerline  and  the  change  in  the 
shapes  of  the  Mach  number  contours  mentioned 
before  is  due  to  these  incoming  and  reflected 
weak  oblique  shock  waves.  Further  evidence  of 
these  waves  can  be  seen  in  Fig.  4,  where  a  plot 
of  the  radial  velocity  along  the  grid  line  j  =  8. 
The  radial  velocity  is  negative  in  the  converg¬ 
ing  section  and  then  it  rapidly  turns  around 
at  the  throat  becoming  positive  in  the  diverg¬ 
ing  section.  The  radial  velocity  then  decreases 
sharply  in  the  neighbourhood  of  z  =  8  and  be¬ 
comes  negative  as  the  flow  traverses  the  inci¬ 
dent  oblique  shock  waves.  The  fluid  is  turned  in 
a  clockwise  direction  towards  the  centerline  by 
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these  waves,  since  they  are  right  running.  The 
radial  velocity  then  increases  just  as  rapidly,  as 
the  flow  traverses  the  oblique  shocks  that  are 
reflected  from  the  centerline.  These  turn  the 
fluid  in  an  anti-clockwise  direction  away  from 
the  centerline,  being  left  running  waves.  The 
radial  velocity  settles  down  to  a  constant  value 
after  about  z  =  14.  The  overturning  of  the 
fluid  as  it  goes  through  the  throat  is  evident 
from  this  figure  when  one  compares  the  value 
of  the  radial  velocity  near  z  «  7  with  the  final 
value  for  z  >  14. 

The  exact  mechanism  by  which  these  com¬ 
pression  waves  are  generated  is  discussed  in  de¬ 
tail  in  [16,  19].  As  the  flow  passes  through  the 
throat,  the  regions  of  the  flow  close  to  the  wall 
are  overturned  so  as  to  be  parallel  to  the  wall. 
Consequently,  when  the  location  where  the  dis¬ 
continuous  change  in  wall  curvature  occurs  is 
reached,  the  flow  must  be  turned  back  towards 
the  centerline,  i.e.  away  from  the  wall.  This 
can  only  occur  through  a  series  of  compression 
waves. 

It  is  likely  that  such  compression  waves  are 
generated  in  the  constrictor  region  of  the  arc- 
jet.  Depending  on  the  physical  shape  of  the 
arc,  it  is  possible  that  flow  turning  such  as 
the  type  just  discussed  occur  at  the  constrictor 
exit.  Another  way  this  effect  can  happen  in  ar- 
ejets  is  via  actual  changes  in  the  contour  of  the 
nozzle  (anode)  wall  due  to  electrode  erosion. 
The  amount  of  erosion  needed  to  trigger  such 
compression  waves  is  rather  small.  Very  sub¬ 
tle  changes  in  the  geometry  can  trigger  these 
compression  waves.  When  a  discharge  is  first 
struck,  the  arc  attaches  either  inside  the  con¬ 
strictor  or  just  at  the  exit  of  the  constrictor. 
As  the  supersonic  flow  becomes  established  and 
steady  state  conditions  are  reached,  the  arc 
attachment  on  the  anode  shifts  in  the  down¬ 
stream  direction.  This  shifting  of  the  anode 
awe  attachment  is  in  reality  accompanied  by 
erosion.  Generation  of  such  compression  waves 
with  possible  coalescence  downstream  into  in¬ 


tersecting  weak  oblique  shocks  has  very  impor¬ 
tant  consequences  for  any  experimental  mea¬ 
surements  that  are  made  at  the  exit  plane. 

4.2  30  kW  Arcjet  Geometry 

Cold  flow  calculations  similar  to  those  de¬ 
scribed  in  the  previous  section,  have  been  per¬ 
formed  for  the  30  KW  arcjet  geometry.  No  ev¬ 
idence  of  compression  waves  is  seen  in  this  ge¬ 
ometry  for  the  case  of  cold  flow.  However,  the 
introduction  of  swirl  at  the  inlet  to  the  thruster 
leads  to  some  interesting  results.  The  profile  of 
the  inlet  swirl  velocity  that  we  use  is  such  that 
it  is  peaked  very  close  to  the  anode  wall,  typ¬ 
ically  within  5  grid  points  (see  Fig.  5).  This 
type  of  swirl  will  be  observed,  for  example,  if 
the  fluid  were  injected  through  ports  periph¬ 
eral  to  the  anode  wall.  Shown  in  Figs.  6  and 
7  are  the  contours  of  swirl  velocities  for  mass 
flow  rates  of  the  order  of  a  few  grams  per  sec 
and  a  few  milligrams  per  second  respectively. 
It  can  be  seen  that  the  introduced  swirl  is  dis¬ 
sipated  in  the  former  while  it  is  enhanced  to 
values  higher  than  the  inlet  values  for  the  latter 
case.  We  find  that  for  the  swirl  to  enhance  in 
the  case  of  higher  mass  flow  rates,  the  plenum 
should  be  shorter. 

The  phenomena  governing  the  enhancement 
or  decay  of  swirl  can  be  explained  from  simple 
considerations  of  conservation  of  angular  mo¬ 
mentum.  Injection  of  swirl  at  the  inlet  repre¬ 
sents  an  influx  of  angular  momentum.  Thus  for 
any  station  downstream,  the  efflux  of  angular 
momentum  will  equal  this  influx  minus  the  to¬ 
tal  torque  due  to  viscous  forces.  If  this  torque 
is  sufficient  to  overcome  the  influx  of  angular 
momentum,  then  the  injected  swirl  decays  in 
the  upstream  plenum  region  itself.  However, 
if  this  torque  is  smaller  than  the  influx  of  an¬ 
gular  momentum,  then  the  swirl  will  decrease 
but  can  persist  until  the  beginning  of  the  con¬ 
verging  section.  Once  this  converging  section  is 
reached,  the  decreasing  area  can  if  the  condi- 
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tions  are  right,  lead  to  an  increase  in  the  swirl 
component  of  the  velocity  from  overall  conser¬ 
vation  of  angular  momentum.  The  total  torque 
is  composed  of  two  quantities.  The  first  is  the 
torque  due  to  the  wall  shear  stresses  (at  anode 
and  cathode),  and  the  second  is  the  torque  due 
to  the  shear  stress  integrated  across  the  flow 
volume  between  cathode  and  anode.  Of  these, 
the  torque  due  to  the  wall  shear  at  the  anode 
is  the  most  important. 

The  critical  parameter  which  determines 
whether  or  not  the  swirl  persists  ans  is  ampli¬ 
fied,  is  the  propellant  mass  flow  rate.  At  higher 
mass  flow  rates  (and  hence  higher  Reynolds 
numbers),  the  boundary  layers  are  closer  to  the 
walls  so  that  the  wall  shear  stresses  are  high. 
Consequently,  the  torque  due  to  viscous  forces 
can  overcome  the  influx  of  angular  momentum 
and  the  swirl  will  perish  long  before  it  reaches 
the  cathode  tip.  In  contrast,  at  the  lower  mass 
flow  rates  (or  lower  Reynolds  numbers),  the 
boundary  layers  are  thicker  so  that  wall  shear 
stresses  are  smaller  leading  to  reduced  torque 
due  to  viscous  forces.  Therefore,  regardless  of 
the  magnitude  of  the  viscosity  there  is  some 
limiting  mass  flow  rate  beyond  which  the  in¬ 
jected  swirl  will  decay.  Obviously,  the  higher 
the  viscosity  the  smaller  the  limiting  mass  flow 
rate.  Thus,  this  mechanism  for  swirl  enhance¬ 
ment  can  indeed  be  operative  at  higher  temper¬ 
atures.  The  mass  flow  rate  and  inlet  plenum 
length  are  thus  key  design  parameters. 

5  Arcjet  Flow 

In  this  section,  we  present  results  for  a  30  kW 
arcjet  geometry,  current  of  50  A,  argon  mass 
flow  rate  of  100  mg/s.  We  wish  to  illustrate 
the  solutions  to  the  one  temperature  model  de¬ 
scribed  in  sections  2  and  3.  We  emphasize  at 
the  outset  that  no  "limits”  or  "floors”  are  ap¬ 
plied  to  the  electrical  conductivity  as  in  [7,  6j. 
The  inlet  stagnation  temperature  is  10000  K 


and  the  inlet  stagnation  pressure  is  10*  Pa.  The 
initial  current  distribution  is  as  shown  in  Fig. 
8.  A  non-dimensional  time  step  of  10-3  (di¬ 
mensional  value  is  5  x  10-9)  is  used  and  conver¬ 
gence  to  atleast  three  decimal  places  is  achieved 
after  2500  time  steps.  The  CPU  time  on  the 
Cray  Y-MP  was  10  minutes. 

Two  cases  were  considered,  one  without  swirl 
and  one  with  the  inlet  swirl  profile  as  shown 
in  Fig.  5.  Figs.  9  through  18,  summarize 
the  results  obtained  after  convergence.  It  can 
be  seen  from  Fig.  9  that  the  current  contours 
have  shifted  downstream  compared  to  the  ini¬ 
tially  prescribed  contours.  This  shift  is  caused 
by  the  varying  conductivity  as  the  supersonic 
flow  becomes  established.  This  reflects  reality 
quite  well  as  when  the  arc  is  first  struck,  it  is 
usuallt  just  downstream  of  the  constrictor  exit 
on  the  anode.  As  the  steady  supersonic  condi¬ 
tions  are  reached,  the  anode  attachment  moves 
downstream.  The  density  of  current  contours 
near  r=0  suggestes  that  the  current  densities 
and  hence  ohmic  heating  are  highest  in  that 
region.  As  expected  the  temperature  and  elec¬ 
tron  number  density  contours  follow  this  very 
well  as  can  be  seen  from  Figs.  11  and  12.  The 
calculated  temperatures  upstream  of  the  cath¬ 
ode  tip  are  in  our  opinion  unrealistically  high. 
This  is  an  inevitable  consequence  of  having  to 
prescribe  a  high  inlet  stagnation  temperature 
in  order  to  enable  the  gas  to  ionize.  This  is 
an  artifact  of  our  simulation  and  caused  by  the 
non- conservative  treatment  of  the  ohmic  heat¬ 
ing  term  in  the  energy  equation.  We  are  at 
present  pursuing  a  conservative  treatment  of 
th  ohmic  heating  term  which  should  yield  much 
lower  temperatures  in  the  inlet  plenum  section. 
The  maximum  temperature  is  33000  K  and  oc¬ 
curs  just  downstream  of  the  cathode  as  can  be 
seen  from  Fig.  11. 

Close  examination  of  the  Mach  contours  in 
Fig.  10  and  the  Mach  profiles  in  Fig.  14,  re¬ 
veal  that  the  flow  reached  supersonic  speeds  in 
the  nozzle  and  becomes  subsonic  by  the  time 
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the  exit  plane  is  reached.  Heat  addition  to  a 
supersonic  flow  is  well  known  to  slow  it  down. 
Examination  of  Figs.  10  and  14  shows  some 
similarity  to  the  cold  flow  solutions  obtained 
for  the  C-D  nozzle  (see  Figs.  2  and  3).  The 
successive  bumps  in  the  Mach  profiles  of  Fig. 
14  suggests  the  generation  of  weak  compression 
waves  which  coalesce  at  the  centerline  form¬ 
ing  an  oblique  shock.  This  complex  series  of 
oblique  shocks  then  reflect  and  exit  the  domain 
via  the  exit  plane.  The  generation  of  such  com¬ 
pression  waves  is  due  to  the  fact  that  once  su¬ 
personic,  the  flow  at  the  exit  of  the  constrictor 
and  beyond  is  forced  to  turn  direction  due  to 
variations  in  ohmic  heating  and  transport  prop¬ 
erties.  Turning  of  supesonic  flows  as  is  well 
known  is  always  accompanied  by  generation  of 
oblique  shocks  or  expansion  fans. 

Contours  of  pressure  are  shown  in  Fig.  13 
and  exhibit  a  decrease  only  until  approximately 
midway  in  the  nozzle.  Although  not  shown 
in  the  figure,  the  pressure  exhibits  an  increase 
near  the  exit  plane  indicating  the  passage  of  the 
flow  through  an  oblique  shock.  Note  the  differ¬ 
ences  in  the  contours  due  different  boundary 
conditions  used  here  vs.  [5].  Fig.  16  shows 
voltage  as  a  function  of  streamwise  distance. 
The  power  expended  in  the  discharge  is  approx¬ 
imately  1  kW.  Fig.  17  shows  contours  of  swirl 
velocity  indicating  that  the  swirl  decays  and 
fails  to  survive  for  this  case.  The  swirl  is  atten¬ 
uated  because  the  torque  due  to  the  shear  stress 
at  the  anode  i£  large.  This  is  an  artifact  of  our 
abnormally  high  inlet  values  of  static  temper¬ 
ature  (approximately  10000  K)  and  ionization 
fractions  (approximately  0.01).  Elimination  of 
the  inlet  plenum  does  not  result  in  amplifica¬ 
tion  of  the  swirl  as  can  be  seen  in  Fig.  18. 
Again  this  is  due  to  the  high  temperatures  and 
ionization  levels  encountered  by  the  flow  before 
the  cathode  tip.  We  must  emphasize  that  we  do 
not  believe  this  to  be  a  true  representation  of 
the  arejet.  Had  we  treated  the  inlet  boundary 
conditions  accounting  for  ohmic  heating,  the 


temperature  and  ionization  level  would  be  far 
lower,  thus  decreasing  the  viscosity  by  orders 
of  magnitude  below  the  level  for  this  case. 

As  far  as  the  effects  of  swirl  are  concerned, 
the  present  cold  flow  results  appear  to  simu¬ 
late  real  arejet  flows  better.  Amplification  of 
swirl  near  the  cathode  tip  would  be  beneficial  in 
providing  much  needed  cooling  to  confine  and 
maintain  the  arc  spot  attachment  area. 

6  Summary  &  Conclusions 

Both  cold  flow  and  ionizing  flow  in  an  arc- 
jet  thruster  have  been  successfully  simulated. 
These  simulations  provide  insight  as  to  the  ef¬ 
fects  of  boundary  conditions  on  the  results  and 
serve  to  identify  the  behavior  of  the  injected 
swirl  at  the  inlet.  It  is  clear  that  swirl  injec¬ 
tion  (both  magnitude  and  location  of  the  max¬ 
imum)  is  intimately  tied  to  two  key  design  pa¬ 
rameters  :  the  operating  mass  flow  rate  and  the 
distance  from  the  inlet  plenum  to  the  beginning 
of  the  converging  section. 

Improvements  are  certainly  in  adequately 
specifying  the  inlet  boundary  conditions,  and 
conservative  treatment  of  the  ohmic  heating 
certainly  would  be  one.  A  two-temperature 
model  along  the  lines  of  Miller  et  al.  [5]  could 
be  an  improvement  for  arejets  operating  on 
monatomic  propellants  or  hydrogen.  However, 
the  two-temperature  treatment  is  suspect  in 
the  case  of  arejets  operating  on  ammonia  or 
hydrazine  due  to  the  presence  of  large  amounts 
of  N2  and  other  diatomic  species  [12].  Fur¬ 
ther  work  incorporating  non- Maxwellian  distri¬ 
butions  of  the  EEDF  would  be  essential  for  am¬ 
monia  and  hydrazine  arejets. 
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30-kW  Aicjet  Geometry 


Schematic  of  the  nominal  30  KW  arcjet  geometry.  The  C-D  nozzle  if  the 
tame  as  this  geometry,  but  without  the  inlet  plenum  or  the  constrictor. 


the  dramatic  flow  turning. 


30  —  IcW  Arcjet  Geometry 
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but  with  no  inlet  plenum. 
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